From 22cfc0bd69c287c404b76b04e536b6c39a3e108e Mon Sep 17 00:00:00 2001
From: Jie Li <76780849+jieli-matrix@users.noreply.github.com>
Date: Fri, 19 Jan 2024 17:38:44 +0800
Subject: [PATCH 01/17] Feature: add python wrapper for math sphbes (#3475)

* recommit for review

* add python wrapper

* remove timer since performace tests add
---
 python/pyabacus/CMakeLists.txt              | 11 +++-
 python/pyabacus/src/py_abacus.cpp           | 13 +++++
 python/pyabacus/src/py_math_base.cpp        | 63 +++++++++++++++++++++
 python/pyabacus/src/py_numerical_radial.cpp |  4 +-
 python/pyabacus/src/pyabacus/__init__.py    |  5 +-
 python/pyabacus/tests/test_base_math.py     | 15 +++++
 python/pyabacus/tests/test_nr.py            | 25 --------
 source/module_base/math_sphbes.cpp          |  7 +--
 8 files changed, 105 insertions(+), 38 deletions(-)
 create mode 100644 python/pyabacus/src/py_abacus.cpp
 create mode 100644 python/pyabacus/src/py_math_base.cpp
 create mode 100644 python/pyabacus/tests/test_base_math.py
 delete mode 100644 python/pyabacus/tests/test_nr.py

diff --git a/python/pyabacus/CMakeLists.txt b/python/pyabacus/CMakeLists.txt
index 399bd4fe57..0effbe83f2 100644
--- a/python/pyabacus/CMakeLists.txt
+++ b/python/pyabacus/CMakeLists.txt
@@ -12,9 +12,14 @@ set(BASE_PATH "${PROJECT_SOURCE_DIR}/../../source/module_base")
 set(ABACUS_SOURCE_DIR "${PROJECT_SOURCE_DIR}/../../source")
 include_directories(${BASE_PATH} ${ABACUS_SOURCE_DIR})
 list(APPEND _sources
-    ${ABACUS_SOURCE_DIR}/module_basis/module_nao/numerical_radial.h
-    ${ABACUS_SOURCE_DIR}/module_basis/module_nao/numerical_radial.cpp
-    ${PROJECT_SOURCE_DIR}/src/py_numerical_radial.cpp)
+    #${ABACUS_SOURCE_DIR}/module_basis/module_nao/numerical_radial.h
+    #${ABACUS_SOURCE_DIR}/module_basis/module_nao/numerical_radial.cpp
+    ${ABACUS_SOURCE_DIR}/module_base/constants.h
+    ${ABACUS_SOURCE_DIR}/module_base/math_sphbes.h
+    ${ABACUS_SOURCE_DIR}/module_base/math_sphbes.cpp
+    ${PROJECT_SOURCE_DIR}/src/py_abacus.cpp
+    #${PROJECT_SOURCE_DIR}/src/py_numerical_radial.cpp
+    ${PROJECT_SOURCE_DIR}/src/py_math_base.cpp)
 python_add_library(_core MODULE ${_sources} WITH_SOABI)
 target_link_libraries(_core PRIVATE pybind11::headers)
 target_compile_definitions(_core PRIVATE VERSION_INFO=${PROJECT_VERSION})
diff --git a/python/pyabacus/src/py_abacus.cpp b/python/pyabacus/src/py_abacus.cpp
new file mode 100644
index 0000000000..34b354dc6b
--- /dev/null
+++ b/python/pyabacus/src/py_abacus.cpp
@@ -0,0 +1,13 @@
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+
+namespace py = pybind11;
+
+void bind_numerical_radial(py::module& m);
+void bind_math_base(py::module& m);
+
+PYBIND11_MODULE(_core, m)
+{
+    // bind_numerical_radial(m);
+    bind_math_base(m);
+}
\ No newline at end of file
diff --git a/python/pyabacus/src/py_math_base.cpp b/python/pyabacus/src/py_math_base.cpp
new file mode 100644
index 0000000000..4378690897
--- /dev/null
+++ b/python/pyabacus/src/py_math_base.cpp
@@ -0,0 +1,63 @@
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+
+#include "module_base/math_sphbes.h"
+
+namespace py = pybind11;
+using namespace pybind11::literals;
+template <typename... Args>
+using overload_cast_ = pybind11::detail::overload_cast_impl<Args...>;
+
+void bind_math_base(py::module& m)
+{
+    py::module module_base = m.def_submodule("ModuleBase");
+
+    py::class_<ModuleBase::Sphbes>(module_base, "Sphbes")
+        .def(py::init<>())
+        .def_static("sphbesj", overload_cast_<const int, const double>()(&ModuleBase::Sphbes::sphbesj), "l"_a, "x"_a)
+        .def_static("dsphbesj", overload_cast_<const int, const double>()(&ModuleBase::Sphbes::dsphbesj), "l"_a, "x"_a)
+        .def_static("sphbesj",
+                    [](const int n, py::array_t<double> r, const double q, const int l, py::array_t<double> jl) {
+                        py::buffer_info r_info = r.request();
+                        if (r_info.ndim != 1)
+                        {
+                            throw std::runtime_error("r array must be 1-dimensional");
+                        }
+                        py::buffer_info jl_info = jl.request();
+                        if (jl_info.ndim != 1)
+                        {
+                            throw std::runtime_error("jl array must be 1-dimensional");
+                        }
+                        ModuleBase::Sphbes::sphbesj(n,
+                                                    static_cast<const double* const>(r_info.ptr),
+                                                    q,
+                                                    l,
+                                                    static_cast<double* const>(jl_info.ptr));
+                    })
+        .def_static("dsphbesj",
+                    [](const int n, py::array_t<double> r, const double q, const int l, py::array_t<double> djl) {
+                        py::buffer_info r_info = r.request();
+                        if (r_info.ndim != 1)
+                        {
+                            throw std::runtime_error("r array must be 1-dimensional");
+                        }
+                        py::buffer_info djl_info = djl.request();
+                        if (djl_info.ndim != 1)
+                        {
+                            throw std::runtime_error("djl array must be 1-dimensional");
+                        }
+                        ModuleBase::Sphbes::dsphbesj(n,
+                                                     static_cast<const double* const>(r_info.ptr),
+                                                     q,
+                                                     l,
+                                                     static_cast<double* const>(djl_info.ptr));
+                    })
+        .def_static("sphbes_zeros", [](const int l, const int n, py::array_t<double> zeros) {
+            py::buffer_info zeros_info = zeros.request();
+            if (zeros_info.ndim != 1)
+            {
+                throw std::runtime_error("zeros array must be 1-dimensional");
+            }
+            ModuleBase::Sphbes::sphbes_zeros(l, n, static_cast<double* const>(zeros_info.ptr));
+        });
+}
\ No newline at end of file
diff --git a/python/pyabacus/src/py_numerical_radial.cpp b/python/pyabacus/src/py_numerical_radial.cpp
index 296229b3d1..ebda8f080b 100644
--- a/python/pyabacus/src/py_numerical_radial.cpp
+++ b/python/pyabacus/src/py_numerical_radial.cpp
@@ -8,7 +8,7 @@ using namespace pybind11::literals;
 template <typename... Args>
 using overload_cast_ = pybind11::detail::overload_cast_impl<Args...>;
 
-PYBIND11_MODULE(_core, m)
+void bind_numerical_radial(py::module& m)
 {
     // Create the submodule for NumericalRadial
     py::module m_numerical_radial = m.def_submodule("NumericalRadial");
@@ -165,4 +165,4 @@ PYBIND11_MODULE(_core, m)
         .def_property_readonly("kgrid", overload_cast_<int>()(&NumericalRadial::kgrid, py::const_))
         .def_property_readonly("rvalue", overload_cast_<int>()(&NumericalRadial::rvalue, py::const_))
         .def_property_readonly("kvalue", overload_cast_<int>()(&NumericalRadial::kvalue, py::const_));
-}
+}
\ No newline at end of file
diff --git a/python/pyabacus/src/pyabacus/__init__.py b/python/pyabacus/src/pyabacus/__init__.py
index cda9318053..94d8c0d5b8 100644
--- a/python/pyabacus/src/pyabacus/__init__.py
+++ b/python/pyabacus/src/pyabacus/__init__.py
@@ -1,3 +1,4 @@
 from __future__ import annotations
-from ._core import __doc__, __version__, NumericalRadial
-__all__ = ["__doc__", "__version__", "NumericalRadial"]
\ No newline at end of file
+# from ._core import __doc__, __version__, NumericalRadial, ModuleBase
+from ._core import ModuleBase
+__all__ = ["ModuleBase"]
\ No newline at end of file
diff --git a/python/pyabacus/tests/test_base_math.py b/python/pyabacus/tests/test_base_math.py
new file mode 100644
index 0000000000..97d5118bac
--- /dev/null
+++ b/python/pyabacus/tests/test_base_math.py
@@ -0,0 +1,15 @@
+from __future__ import annotations
+
+import pyabacus as m
+import numpy as np
+
+
+def test_version():
+    assert m.__version__ == "0.0.1"
+
+def test_sphbes():
+    s = m.ModuleBase.Sphbes()
+    # test for sphbesj
+    assert s.sphbesj(1, 0.0) == 0.0
+    assert s.sphbesj(0, 0.0) == 1.0
+
diff --git a/python/pyabacus/tests/test_nr.py b/python/pyabacus/tests/test_nr.py
deleted file mode 100644
index 4986331b25..0000000000
--- a/python/pyabacus/tests/test_nr.py
+++ /dev/null
@@ -1,25 +0,0 @@
-from __future__ import annotations
-
-import pyabacus as m
-
-
-def test_version():
-    assert m.__version__ == "0.0.1"
-
-def test_attributes():
-    chi = m.NumericalRadial()
-    # string
-    assert chi.symbol == ''
-    # integer
-    assert chi.itype == 0
-    assert chi.izeta == 0
-    assert chi.l == -1
-    assert chi.nr == 0
-    assert chi.nk == 0
-    # float
-    assert chi.rcut == 0.0
-    assert chi.kcut == 0.0
-    assert chi.pr == 0.0
-    assert chi.pk == 0.0
-    # bool
-    assert chi.is_fft_compliant == False
diff --git a/source/module_base/math_sphbes.cpp b/source/module_base/math_sphbes.cpp
index 00326f928e..73e0127e6b 100644
--- a/source/module_base/math_sphbes.cpp
+++ b/source/module_base/math_sphbes.cpp
@@ -1,7 +1,7 @@
 #include "math_sphbes.h"
-#include "timer.h"
 #include "constants.h"
 #include <algorithm>
+#include <iostream>
 
 #include <cassert>
 
@@ -425,7 +425,6 @@ void Sphbes::Spherical_Bessel
     double *jl		 // jl(1:msh) = j_l(q*r(i)),spherical bessel function
 )
 {
-    ModuleBase::timer::tick("Sphbes","Spherical_Bessel");
     double x1=0.0;
 
     int i=0;
@@ -598,7 +597,6 @@ void Sphbes::Spherical_Bessel
         }
     }
 
-    ModuleBase::timer::tick("Sphbes","Spherical_Bessel");
     return;
 }
 
@@ -613,7 +611,6 @@ void Sphbes::Spherical_Bessel
 	double *sjp
 )
 {
-	ModuleBase::timer::tick("Sphbes","Spherical_Bessel");
 
 	//calculate jlx first
 	Spherical_Bessel (msh, r, q, l, sj);
@@ -634,7 +631,6 @@ void Sphbes::dSpherical_Bessel_dx
     double *djl		 // jl(1:msh) = j_l(q*r(i)),spherical bessel function
 )
 {
-    ModuleBase::timer::tick("Sphbes","dSpherical_Bessel_dq");
     if (l < 0 )
     {
 		std::cout << "We temporarily only calculate derivative of l >= 0." << std::endl;
@@ -682,7 +678,6 @@ void Sphbes::dSpherical_Bessel_dx
         }
         delete[] jl;
     }
-    ModuleBase::timer::tick("Sphbes","dSpherical_Bessel_dq");
     return;
 }
 

From 115a9899a6a95fcc3fc97566911c3ad9c014048a Mon Sep 17 00:00:00 2001
From: kirk0830 <67682086+kirk0830@users.noreply.github.com>
Date: Fri, 19 Jan 2024 18:00:54 +0800
Subject: [PATCH 02/17] Feature: support segment split in kline mode in KPT
 file and `out_band` band output precision control, `8` as default (#3493)

* add precision control

* correct serial version of nscf_band function

* fix issue 3482

* update unit and integrated test

* update document

* correct unittest and make compatible with false and true
---
 docs/advanced/input_files/input-main.md       |  4 +-
 source/module_cell/klist.cpp                  | 44 +++++++++++++------
 source/module_cell/klist.h                    |  1 +
 source/module_esolver/esolver_ks_lcao.cpp     | 14 ++++--
 source/module_esolver/esolver_ks_pw.cpp       |  5 ++-
 source/module_io/input.cpp                    | 13 +++---
 source/module_io/input.h                      | 12 ++++-
 source/module_io/nscf_band.cpp                | 41 ++++++++++++-----
 source/module_io/nscf_band.h                  |  1 +
 source/module_io/parameter_pool.cpp           |  2 +-
 source/module_io/test/input_test.cpp          | 14 ++++--
 source/module_io/test/input_test_para.cpp     |  3 +-
 source/module_io/test/support/INPUT           |  2 +-
 source/module_io/test/support/witestfile      |  2 +-
 .../module_io/test_serial/nscf_band_test.cpp  | 11 ++++-
 source/module_io/write_input.cpp              |  4 +-
 .../107_PW_OBOD_MemSaver/refBANDS_1.dat       | 12 ++---
 .../107_PW_OB_outputbands/refBANDS_1.dat      | 12 ++---
 tests/integrate/207_NO_KP_OB/refBANDS_1.dat   | 12 ++---
 19 files changed, 140 insertions(+), 69 deletions(-)

diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
index b9ee122af6..fb721a7486 100644
--- a/docs/advanced/input_files/input-main.md
+++ b/docs/advanced/input_files/input-main.md
@@ -1494,8 +1494,8 @@ These variables are used to control the output of properties.
 
 ### out_band
 
-- **Type**: Boolean
-- **Description**: Whether to output the band structure (in eV). For more information, refer to the [band.md](../elec_properties/band.md)
+- **Type**: Boolean Integer(optional)
+- **Description**: Whether to output the band structure (in eV), optionally output precision can be set by a second parameter, default is 8. For more information, refer to the [band.md](../elec_properties/band.md)
 - **Default**: False
 
 ### out_proj_band
diff --git a/source/module_cell/klist.cpp b/source/module_cell/klist.cpp
index 52bc42440d..e11ab183f5 100644
--- a/source/module_cell/klist.cpp
+++ b/source/module_cell/klist.cpp
@@ -362,6 +362,10 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 
 			//recalculate nkstot.
 			nkstot = 0;
+            /* ISSUE#3482: to distinguish different kline segments */
+            std::vector<int> kpt_segids;
+            kl_segids.clear(); kl_segids.shrink_to_fit();
+            int kpt_segid = 0;
 			for(int iks=0; iks<nks_special; iks++)
 			{
 				ifk >> ksx[iks];
@@ -371,6 +375,9 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 				//std::cout << " nkl[" << iks << "]=" << nkl[iks] << std::endl;
 				assert(nkl[iks] >= 0);
 				nkstot += nkl[iks];
+                /* ISSUE#3482: to distinguish different kline segments */
+                if((nkl[iks] == 1)&&(iks!=(nks_special-1))) kpt_segid++;
+                kpt_segids.push_back(kpt_segid);
 			}
 			assert( nkl[nks_special-1] == 1);
 
@@ -389,6 +396,7 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 					kvec_c[count].x = ksx[iks-1] + is*dx;
 					kvec_c[count].y = ksy[iks-1] + is*dy;
 					kvec_c[count].z = ksz[iks-1] + is*dz;
+                    kl_segids.push_back(kpt_segids[iks-1]); /* ISSUE#3482: to distinguish different kline segments */
 					++count;
 				}
 			}
@@ -397,15 +405,14 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 			kvec_c[count].x = ksx[nks_special-1];
 			kvec_c[count].y = ksy[nks_special-1];
 			kvec_c[count].z = ksz[nks_special-1];
+            kl_segids.push_back(kpt_segids[nks_special-1]); /* ISSUE#3482: to distinguish different kline segments */
 			++count;
 
 			//std::cout << " count = " << count << std::endl;
-			assert (count == nkstot );
-
-			for(int ik=0; ik<nkstot; ik++)
-			{
-				wk[ik] = 1.0;
-			}
+			assert(count == nkstot);
+            assert(kl_segids.size() == nkstot); /* ISSUE#3482: to distinguish different kline segments */
+			
+            std::for_each(wk.begin(), wk.end(), [](double& d){d = 1.0;});
 
             this->kc_done = true;
 
@@ -439,15 +446,22 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 
 			//recalculate nkstot.
 			nkstot = 0;
+            /* ISSUE#3482: to distinguish different kline segments */
+            std::vector<int> kpt_segids;
+            kl_segids.clear(); kl_segids.shrink_to_fit();
+            int kpt_segid = 0;
 			for(int iks=0; iks<nks_special; iks++)
 			{
 				ifk >> ksx[iks];
 				ifk >> ksy[iks];
 				ifk >> ksz[iks];
-				ModuleBase::GlobalFunc::READ_VALUE( ifk, nkl[iks] );
+				ModuleBase::GlobalFunc::READ_VALUE( ifk, nkl[iks] ); /* so ifk is ifstream for kpoint, then nkl is number of kpoints on line */
 				//std::cout << " nkl[" << iks << "]=" << nkl[iks] << std::endl;
 				assert(nkl[iks] >= 0);
 				nkstot += nkl[iks];
+                /* ISSUE#3482: to distinguish different kline segments */
+                if((nkl[iks] == 1)&&(iks!=(nks_special-1))) kpt_segid++;
+                kpt_segids.push_back(kpt_segid);
 			}
 			assert( nkl[nks_special-1] == 1);
 
@@ -466,6 +480,7 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 					kvec_d[count].x = ksx[iks-1] + is*dx;
 					kvec_d[count].y = ksy[iks-1] + is*dy;
 					kvec_d[count].z = ksz[iks-1] + is*dz;
+                    kl_segids.push_back(kpt_segids[iks-1]); /* ISSUE#3482: to distinguish different kline segments */
 					++count;
 				}
 			}
@@ -474,18 +489,16 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 			kvec_d[count].x = ksx[nks_special-1];
 			kvec_d[count].y = ksy[nks_special-1];
 			kvec_d[count].z = ksz[nks_special-1];
+            kl_segids.push_back(kpt_segids[nks_special-1]); /* ISSUE#3482: to distinguish different kline segments */
 			++count;
 
 			//std::cout << " count = " << count << std::endl;
-			assert (count == nkstot );
+			assert(count == nkstot );
+            assert(kl_segids.size() == nkstot); /* ISSUE#3482: to distinguish different kline segments */
 
-			for(int ik=0; ik<nkstot; ik++)
-			{
-				wk[ik] = 1.0;
-			}
+			std::for_each(wk.begin(), wk.end(), [](double& d){d = 1.0;});
 
             this->kd_done = true;
-
 		}
 
         else
@@ -1122,6 +1135,9 @@ void K_Vectors::mpi_k(void)
 
     Parallel_Common::bcast_int(nmp, 3);
 
+    kl_segids.resize(nkstot);
+    Parallel_Common::bcast_int(kl_segids.data(), nkstot);
+
     Parallel_Common::bcast_double(koffset, 3);
 
     this->nks = GlobalC::Pkpoints.nks_pool[GlobalV::MY_POOL];
@@ -1352,6 +1368,8 @@ void K_Vectors::mpi_k_after_vc(void)
     Parallel_Common::bcast_int(nspin);
     Parallel_Common::bcast_int(nkstot);
     Parallel_Common::bcast_int(nmp, 3);
+    kl_segids.resize(nkstot);
+    Parallel_Common::bcast_int(kl_segids.data(), nkstot);
     Parallel_Common::bcast_double(koffset, 3);
 
     this->nks = GlobalC::Pkpoints.nks_pool[GlobalV::MY_POOL];
diff --git a/source/module_cell/klist.h b/source/module_cell/klist.h
index a9e06f8614..aa92cf29fd 100644
--- a/source/module_cell/klist.h
+++ b/source/module_cell/klist.h
@@ -29,6 +29,7 @@ class K_Vectors
     int nkstot_full;    /// number of k points in full k mesh
 
     int nmp[3];						// Number of Monhorst-Pack
+    std::vector<int> kl_segids;	// index of kline segment
 
     K_Vectors();
     ~K_Vectors();
diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
index 8fdc9217fd..a7d0a5abdd 100644
--- a/source/module_esolver/esolver_ks_lcao.cpp
+++ b/source/module_esolver/esolver_ks_lcao.cpp
@@ -309,7 +309,7 @@ namespace ModuleESolver
     GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl;
     GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl;
 
-    if (INPUT.out_dos != 0 || INPUT.out_band != 0 || INPUT.out_proj_band != 0)
+    if (INPUT.out_dos != 0 || INPUT.out_band[0] != 0 || INPUT.out_proj_band != 0)
     {
         GlobalV::ofs_running << "\n\n\n\n";
         GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
@@ -331,7 +331,7 @@ namespace ModuleESolver
 
     int nspin0 = (GlobalV::NSPIN == 2) ? 2 : 1;
 
-    if (INPUT.out_band) // pengfei 2014-10-13
+    if (INPUT.out_band[0]) // pengfei 2014-10-13
     {
         int nks = 0;
         if (nspin0 == 1)
@@ -348,7 +348,15 @@ namespace ModuleESolver
             std::stringstream ss2;
             ss2 << GlobalV::global_out_dir << "BANDS_" << is + 1 << ".dat";
             GlobalV::ofs_running << "\n Output bands in file: " << ss2.str() << std::endl;
-            ModuleIO::nscf_band(is, ss2.str(), nks, GlobalV::NBANDS, 0.0, this->pelec->ekb, this->kv, &(GlobalC::Pkpoints));
+            ModuleIO::nscf_band(is, 
+                                ss2.str(), 
+                                nks, 
+                                GlobalV::NBANDS, 
+                                0.0, 
+                                INPUT.out_band[1],
+                                this->pelec->ekb, 
+                                this->kv, 
+                                &(GlobalC::Pkpoints));
         }
     } // out_band
 
diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
index 294f1636f7..2e81938ae2 100644
--- a/source/module_esolver/esolver_ks_pw.cpp
+++ b/source/module_esolver/esolver_ks_pw.cpp
@@ -959,7 +959,7 @@ void ESolver_KS_PW<T, Device>::postprocess()
     GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl;
     GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl;
 
-    if (INPUT.out_dos != 0 || INPUT.out_band != 0)
+    if (INPUT.out_dos != 0 || INPUT.out_band[0] != 0)
     {
         GlobalV::ofs_running << "\n\n\n\n";
         GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
@@ -1001,7 +1001,7 @@ void ESolver_KS_PW<T, Device>::postprocess()
         }
     }
 
-    if (INPUT.out_band) // pengfei 2014-10-13
+    if (INPUT.out_band[0]) // pengfei 2014-10-13
     {
         int nks = 0;
         if (nspin0 == 1)
@@ -1022,6 +1022,7 @@ void ESolver_KS_PW<T, Device>::postprocess()
                                 nks,
                                 GlobalV::NBANDS,
                                 0.0,
+                                INPUT.out_band[1],
                                 this->pelec->ekb,
                                 this->kv,
                                 &(GlobalC::Pkpoints));
diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp
index 2685ba540e..abe8f28ebd 100644
--- a/source/module_io/input.cpp
+++ b/source/module_io/input.cpp
@@ -336,7 +336,7 @@ void Input::Default(void)
     out_wfc_pw = 0;
     out_wfc_r = 0;
     out_dos = 0;
-    out_band = 0;
+    out_band = {0, 8};
     out_proj_band = 0;
     out_mat_hs = {0, 8};
     out_mat_xc = 0;
@@ -1378,13 +1378,13 @@ bool Input::Read(const std::string& fn)
         }
         else if (strcmp("out_band", word) == 0)
         {
-            read_bool(ifs, out_band);
+            read_value2stdvector(ifs, out_band);
+            if(out_band.size() == 1) out_band.push_back(8);
         }
         else if (strcmp("out_proj_band", word) == 0)
         {
             read_bool(ifs, out_proj_band);
         }
-
         else if (strcmp("out_mat_hs", word) == 0)
         {
             read_value2stdvector(ifs, out_mat_hs);
@@ -2826,7 +2826,7 @@ void Input::Default_2(void) // jiyy add 2019-08-04
         this->relax_nmax = 1;
         out_stru = 0;
         out_dos = 0;
-        out_band = 0;
+        out_band[0] = 0;
         out_proj_band = 0;
         cal_force = 0;
         init_wfc = "file";
@@ -2843,7 +2843,7 @@ void Input::Default_2(void) // jiyy add 2019-08-04
         this->relax_nmax = 1;
         out_stru = 0;
         out_dos = 0;
-        out_band = 0;
+        out_band[0] = 0;
         out_proj_band = 0;
         cal_force = 0;
         init_wfc = "file";
@@ -3325,7 +3325,8 @@ void Input::Bcast()
     Parallel_Common::bcast_int(out_wfc_pw);
     Parallel_Common::bcast_bool(out_wfc_r);
     Parallel_Common::bcast_int(out_dos);
-    Parallel_Common::bcast_bool(out_band);
+    if(GlobalV::MY_RANK != 0) out_band.resize(2); /* If this line is absent, will cause segmentation fault in io_input_test_para */
+    Parallel_Common::bcast_int(out_band.data(), 2);
     Parallel_Common::bcast_bool(out_proj_band);
     if(GlobalV::MY_RANK != 0) out_mat_hs.resize(2); /* If this line is absent, will cause segmentation fault in io_input_test_para */
     Parallel_Common::bcast_int(out_mat_hs.data(), 2);
diff --git a/source/module_io/input.h b/source/module_io/input.h
index b4e983abad..6393483cb4 100644
--- a/source/module_io/input.h
+++ b/source/module_io/input.h
@@ -263,7 +263,7 @@ class Input
     int out_wfc_pw; // 0: no; 1: txt; 2: dat
     bool out_wfc_r; // 0: no; 1: yes
     int out_dos; // dos calculation. mohan add 20090909
-    bool out_band; // band calculation pengfei 2014-10-13
+    std::vector<int> out_band; // band calculation pengfei 2014-10-13
     bool out_proj_band; // projected band structure calculation jiyy add 2022-05-11
     std::vector<int> out_mat_hs; // output H matrix and S matrix in local basis.
     bool out_mat_xc; // output exchange-correlation matrix in KS-orbital representation.
@@ -667,7 +667,15 @@ class Input
     template <typename T>
     typename std::enable_if<std::is_same<T, double>::value, T>::type cast_string(const std::string& str) { return std::stod(str); }
     template <typename T>
-    typename std::enable_if<std::is_same<T, int>::value, T>::type cast_string(const std::string& str) { return std::stoi(str); }
+    typename std::enable_if<std::is_same<T, int>::value, T>::type cast_string(const std::string& str)
+    {
+        if (str == "true" || str == "1")
+            return 1;
+        else if (str == "false" || str == "0")
+            return 0;
+        else
+            return std::stoi(str);
+    }
     template <typename T>
     typename std::enable_if<std::is_same<T, bool>::value, T>::type cast_string(const std::string& str) { return (str == "true" || str == "1"); }
     template <typename T>
diff --git a/source/module_io/nscf_band.cpp b/source/module_io/nscf_band.cpp
index d8b7b05ca6..290dc58bd3 100644
--- a/source/module_io/nscf_band.cpp
+++ b/source/module_io/nscf_band.cpp
@@ -3,6 +3,7 @@
 #include "module_base/global_variable.h"
 #include "module_base/timer.h"
 #include "module_base/tool_title.h"
+#include "module_base/formatter_physfmt.h"
 
 void ModuleIO::nscf_band(
 	const int &is,
@@ -10,6 +11,7 @@ void ModuleIO::nscf_band(
 	const int &nks, 
 	const int &nband,
 	const double &fermie,
+	const int &precision,
 	const ModuleBase::matrix& ekb,
 	const K_Vectors& kv,
 	const Parallel_Kpoints* Pkpoints)
@@ -33,23 +35,28 @@ void ModuleIO::nscf_band(
 		if (ik>0)
 		{
 			auto delta=kv.kvec_c[ik]-kv.kvec_c[ik-1];
-			klength[ik] = klength[ik-1] + delta.norm();
+			klength[ik] = klength[ik-1];
+			klength[ik] += (kv.kl_segids[ik] == kv.kl_segids[ik-1]) ? delta.norm() : 0.0;
 		}
+		/* first find if present kpoint in present pool */
 		if ( GlobalV::MY_POOL == Pkpoints->whichpool[ik] )
 		{
+			/* then get the local kpoint index, which starts definitly from 0 */
 			const int ik_now = ik - Pkpoints->startk_pool[GlobalV::MY_POOL];
+			/* if present kpoint corresponds the spin of the present one */
 			if( kv.isk[ik_now+is*nks] == is )
 			{ 
 				if ( GlobalV::RANK_IN_POOL == 0)
 				{
-					std::ofstream ofs(out_band_dir.c_str(),std::ios::app);
-					ofs << std::setprecision(8);
-					//start from 1
-					ofs << ik+1;
-					ofs << " " << klength[ik] << " ";
+					formatter::PhysicalFmt physfmt; // create a physical formatter temporarily
+					std::ofstream ofs(out_band_dir.c_str(), std::ios::app);
+					physfmt.adjust_formatter_flexible(4, 0, false); // for integer
+					ofs << physfmt.get_p_formatter()->format(ik+1);
+					physfmt.adjust_formatter_flexible(precision, 4.0/double(precision), false); // for decimal
+					ofs << physfmt.get_p_formatter()->format(klength[ik]);
 					for(int ib = 0; ib < nband; ib++)
 					{
-						ofs << " " << (ekb(ik_now+is*nks, ib)-fermie) * ModuleBase::Ry_to_eV;
+						ofs << physfmt.get_p_formatter()->format((ekb(ik_now+is*nks, ib)-fermie) * ModuleBase::Ry_to_eV);
 					}
 					ofs << std::endl;
 					ofs.close();	
@@ -83,18 +90,30 @@ void ModuleIO::nscf_band(
 #else
 //	std::cout<<"\n nband = "<<nband<<std::endl;
 //	std::cout<<out_band_dir<<std::endl;
-
+	formatter::PhysicalFmt physfmt; // create a physical formatter temporarily
+	std::vector<double> klength;
+	klength.resize(nks);
+	klength[0] = 0.0;
 	std::ofstream ofs(out_band_dir.c_str());
 	for(int ik=0;ik<nks;ik++)
 	{
+		if (ik>0)
+		{
+			auto delta=kv.kvec_c[ik]-kv.kvec_c[ik-1];
+			klength[ik] = klength[ik-1];
+			klength[ik] += (kv.kl_segids[ik] == kv.kl_segids[ik-1]) ? delta.norm() : 0.0;
+		}
 		if( kv.isk[ik] == is)
 		{
-			ofs<<std::setw(12)<<ik + 1;
+			physfmt.adjust_formatter_flexible(4, 0, false); // for integer
+			ofs << physfmt.get_p_formatter()->format(ik+1);
+			physfmt.adjust_formatter_flexible(precision, 4.0/double(precision), false); // for decimal
+			ofs << physfmt.get_p_formatter()->format(klength[ik]); // add klength, in accordance with the MPI version
 			for(int ibnd = 0; ibnd < nband; ibnd++)
 			{
-				ofs <<std::setw(15) << (ekb(ik, ibnd)-fermie) * ModuleBase::Ry_to_eV;
+				ofs << physfmt.get_p_formatter()->format((ekb(ik, ibnd)-fermie) * ModuleBase::Ry_to_eV);
 			}
-			ofs<<std::endl;
+			ofs << std::endl;
 		}
 	}
 	ofs.close();
diff --git a/source/module_io/nscf_band.h b/source/module_io/nscf_band.h
index 6a22427551..3ec96d4a9f 100644
--- a/source/module_io/nscf_band.h
+++ b/source/module_io/nscf_band.h
@@ -12,6 +12,7 @@ namespace ModuleIO
 		const int &nks, 
 		const int &nband, 
 		const double &fermie,
+		const int &precision,
 		const ModuleBase::matrix &ekb,
 		const K_Vectors& kv,
 		const Parallel_Kpoints* Pkpoints);
diff --git a/source/module_io/parameter_pool.cpp b/source/module_io/parameter_pool.cpp
index 906b9b57d8..9a6fbed3f4 100644
--- a/source/module_io/parameter_pool.cpp
+++ b/source/module_io/parameter_pool.cpp
@@ -917,7 +917,7 @@ bool input_parameters_set(std::map<std::string, InputParameter> input_parameters
     }
     else if (input_parameters.count("out_band") != 0)
     {
-        INPUT.out_band = *static_cast<bool*>(input_parameters["out_band"].get());
+        INPUT.out_band = *static_cast<std::vector<int>*>(input_parameters["out_band"].get());
     }
     else if (input_parameters.count("out_proj_band") != 0)
     {
diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp
index 02a5a19e10..11bce873ab 100644
--- a/source/module_io/test/input_test.cpp
+++ b/source/module_io/test/input_test.cpp
@@ -176,9 +176,11 @@ TEST_F(InputTest, Default)
         EXPECT_EQ(INPUT.out_wfc_pw,0);
         EXPECT_EQ(INPUT.out_wfc_r,0);
         EXPECT_EQ(INPUT.out_dos,0);
-        EXPECT_EQ(INPUT.out_band,0);
+        EXPECT_EQ(INPUT.out_band[0],0);
+		EXPECT_EQ(INPUT.out_band[1],8);
         EXPECT_EQ(INPUT.out_proj_band,0);
         EXPECT_EQ(INPUT.out_mat_hs[0],0);
+		EXPECT_EQ(INPUT.out_mat_hs[1],8);
         EXPECT_EQ(INPUT.out_mat_hs2,0);
         EXPECT_EQ(INPUT.out_mat_xc, 0);
         EXPECT_EQ(INPUT.out_interval,1);
@@ -539,9 +541,11 @@ TEST_F(InputTest, Read)
         EXPECT_EQ(INPUT.out_wfc_pw,0);
         EXPECT_EQ(INPUT.out_wfc_r,0);
         EXPECT_EQ(INPUT.out_dos,0);
-        EXPECT_EQ(INPUT.out_band,0);
+        EXPECT_EQ(INPUT.out_band[0],0);
+		EXPECT_EQ(INPUT.out_band[1],8);
         EXPECT_EQ(INPUT.out_proj_band,0);
         EXPECT_EQ(INPUT.out_mat_hs[0],0);
+		EXPECT_EQ(INPUT.out_mat_hs[1],8);
         EXPECT_EQ(INPUT.out_mat_hs2,0);
         EXPECT_EQ(INPUT.out_mat_xc, 0);
         EXPECT_EQ(INPUT.out_interval,1);
@@ -921,7 +925,8 @@ TEST_F(InputTest, Default_2)
     EXPECT_EQ(INPUT.relax_nmax, 1);
     EXPECT_EQ(INPUT.out_stru, 0);
     EXPECT_EQ(INPUT.symmetry, "0");
-	EXPECT_EQ(INPUT.out_band,0);
+	EXPECT_EQ(INPUT.out_band[0],0);
+	EXPECT_EQ(INPUT.out_band[1],8);
 	EXPECT_EQ(INPUT.out_proj_band,0);
 	EXPECT_EQ(INPUT.cal_force,0);
 	EXPECT_EQ(INPUT.init_wfc,"file");
@@ -943,7 +948,8 @@ TEST_F(InputTest, Default_2)
     EXPECT_EQ(INPUT.relax_nmax, 1);
     EXPECT_EQ(INPUT.symmetry, "0");
     EXPECT_EQ(INPUT.out_stru, 0);
-	EXPECT_EQ(INPUT.out_band,0);
+	EXPECT_EQ(INPUT.out_band[0],0);
+	EXPECT_EQ(INPUT.out_band[1],8);
 	EXPECT_EQ(INPUT.out_proj_band,0);
 	EXPECT_EQ(INPUT.cal_force,0);
 	EXPECT_EQ(INPUT.init_wfc,"file");
diff --git a/source/module_io/test/input_test_para.cpp b/source/module_io/test/input_test_para.cpp
index 8ae4bf63d1..7cc4d6bc73 100644
--- a/source/module_io/test/input_test_para.cpp
+++ b/source/module_io/test/input_test_para.cpp
@@ -181,7 +181,8 @@ TEST_F(InputParaTest, Bcast)
     EXPECT_EQ(INPUT.out_wfc_pw, 0);
     EXPECT_EQ(INPUT.out_wfc_r, 0);
     EXPECT_EQ(INPUT.out_dos, 0);
-    EXPECT_EQ(INPUT.out_band, 0);
+    EXPECT_EQ(INPUT.out_band[0], 0);
+    EXPECT_EQ(INPUT.out_band[1], 8);
     EXPECT_EQ(INPUT.out_proj_band, 0);
     EXPECT_EQ(INPUT.out_mat_hs[0], 0);
     EXPECT_EQ(INPUT.out_mat_hs[1], 8);
diff --git a/source/module_io/test/support/INPUT b/source/module_io/test/support/INPUT
index 469dff2ff4..4fbde867db 100644
--- a/source/module_io/test/support/INPUT
+++ b/source/module_io/test/support/INPUT
@@ -59,7 +59,7 @@ out_pot                        2 #output realspace potential
 out_wfc_pw                     0 #output wave functions
 out_wfc_r                      0 #output wave functions in realspace
 out_dos                        0 #output energy and dos
-out_band                       false #output energy and band structure
+out_band                       0 #output energy and band structure
 out_proj_band                  FaLse #output projected band structure
 restart_save                   f #print to disk every step for restart
 restart_load                   F #restart from disk
diff --git a/source/module_io/test/support/witestfile b/source/module_io/test/support/witestfile
index 4043773876..4db819d53f 100644
--- a/source/module_io/test/support/witestfile
+++ b/source/module_io/test/support/witestfile
@@ -55,7 +55,7 @@ out_pot                        2 #output realspace potential
 out_wfc_pw                     0 #output wave functions
 out_wfc_r                      0 #output wave functions in realspace
 out_dos                        0 #output energy and dos
-out_band                       false #output energy and band structure
+out_band                       0 #output energy and band structure
 out_proj_band                  FaLse #output projected band structure
 restart_save                   f #print to disk every step for restart
 restart_load                   F #restart from disk
diff --git a/source/module_io/test_serial/nscf_band_test.cpp b/source/module_io/test_serial/nscf_band_test.cpp
index 4483bf37cd..db9bf752fb 100644
--- a/source/module_io/test_serial/nscf_band_test.cpp
+++ b/source/module_io/test_serial/nscf_band_test.cpp
@@ -54,9 +54,16 @@ class BandTest : public ::testing::Test
 	    ekb(1,1) =  2.0;
 	    ekb(1,2) =  3.0;
         kv = new K_Vectors;
+        // specify the kpoints
+        kv->kvec_c.resize(nks);
+        kv->kvec_c[0] = ModuleBase::Vector3<double>(0.0, 0.0, 0.0);
+        kv->kvec_c[1] = ModuleBase::Vector3<double>(1.0, 0.0, 0.0);
         kv->isk.resize(nks);
         kv->isk[0] = 0;
         kv->isk[1] = 1;
+        kv->kl_segids.resize(nks);
+        kv->kl_segids[0] = 0;
+        kv->kl_segids[1] = 0;
         Pkpoints = new Parallel_Kpoints;
     }
 
@@ -81,12 +88,12 @@ class BandTest : public ::testing::Test
 TEST_F(BandTest, nscf_band)
 {
     // Call the function to be tested
-    ModuleIO::nscf_band(is, out_band_dir, nks, nband, fermie, ekb, *kv, Pkpoints);
+    ModuleIO::nscf_band(is, out_band_dir, nks, nband, fermie, 8, ekb, *kv, Pkpoints);
 
     // Check the output file
     std::ifstream ifs(out_band_dir);
     std::string str((std::istreambuf_iterator<char>(ifs)),std::istreambuf_iterator<char>());
     ASSERT_TRUE(ifs.is_open());
-    EXPECT_THAT(str, testing::HasSubstr("1       -27.2114       -13.6057              0"));
+    EXPECT_THAT(str, testing::HasSubstr("1   0.00000000 -27.21139600 -13.60569800   0.00000000"));
     ifs.close();
 }
diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp
index 5c0fc4cada..b7e49901c8 100644
--- a/source/module_io/write_input.cpp
+++ b/source/module_io/write_input.cpp
@@ -122,7 +122,7 @@ void Input::Print(const std::string &fn) const
     ModuleBase::GlobalFunc::OUTP(ofs, "out_wfc_pw", out_wfc_pw, "output wave functions");
     ModuleBase::GlobalFunc::OUTP(ofs, "out_wfc_r", out_wfc_r, "output wave functions in realspace");
     ModuleBase::GlobalFunc::OUTP(ofs, "out_dos", out_dos, "output energy and dos");
-    ModuleBase::GlobalFunc::OUTP(ofs, "out_band", out_band, "output energy and band structure");
+    ModuleBase::GlobalFunc::OUTP(ofs, "out_band", out_band[0], "output energy and band structure (with precision "+std::to_string(out_band[1])+")");
     ModuleBase::GlobalFunc::OUTP(ofs, "out_proj_band", out_proj_band, "output projected band structure");
     ModuleBase::GlobalFunc::OUTP(ofs, "restart_save", restart_save, "print to disk every step for restart");
     ModuleBase::GlobalFunc::OUTP(ofs, "restart_load", restart_load, "restart from disk");
@@ -222,7 +222,7 @@ ModuleBase::GlobalFunc::OUTP(ofs, "out_bandgap", out_bandgap, "if true, print ou
     ModuleBase::GlobalFunc::OUTP(ofs, "lcao_dk", lcao_dk, "delta k for 1D integration in LCAO");
     ModuleBase::GlobalFunc::OUTP(ofs, "lcao_dr", lcao_dr, "delta r for 1D integration in LCAO");
     ModuleBase::GlobalFunc::OUTP(ofs, "lcao_rmax", lcao_rmax, "max R for 1D two-center integration table");
-    ModuleBase::GlobalFunc::OUTP(ofs, "out_mat_hs", out_mat_hs[0], "output H and S matrix");
+    ModuleBase::GlobalFunc::OUTP(ofs, "out_mat_hs", out_mat_hs[0], "output H and S matrix (with precision "+std::to_string(out_mat_hs[1])+")");
     ModuleBase::GlobalFunc::OUTP(ofs, "out_mat_hs2", out_mat_hs2, "output H(R) and S(R) matrix");
     ModuleBase::GlobalFunc::OUTP(ofs, "out_mat_dh", out_mat_dh, "output of derivative of H(R) matrix");
     ModuleBase::GlobalFunc::OUTP(ofs, "out_mat_xc", out_mat_xc, "output exchange-correlation matrix in KS-orbital representation");
diff --git a/tests/integrate/107_PW_OBOD_MemSaver/refBANDS_1.dat b/tests/integrate/107_PW_OBOD_MemSaver/refBANDS_1.dat
index af0ad58c0d..1aa4b94ac1 100644
--- a/tests/integrate/107_PW_OBOD_MemSaver/refBANDS_1.dat
+++ b/tests/integrate/107_PW_OBOD_MemSaver/refBANDS_1.dat
@@ -1,6 +1,6 @@
-1 0  -3.3870593 -0.79801307 5.0648821 5.0648821 7.8411435 9.605949
-2 0.17320508  -3.8620194 -0.075179882 5.1365314 5.1365314 7.9181006 9.6849562
-3 0.34641016  -4.6307963 1.434193 5.3528485 5.3528485 8.1554014 9.8149803
-4 0.51961524  -5.2581292 3.25095 5.6954121 5.6954121 8.5186442 9.653424
-5 0.69282032  -5.6519805 5.1370698 6.0846605 6.0846605 8.8678468 9.1370653
-6 0.8660254  -5.7858738 6.2887673 6.2887673 6.2887673 8.8364355 8.8364355
+   1   0.00000000  -3.38705933  -0.79801307   5.06488210   5.06488210   7.84114355   9.60594903
+   2   0.17320508  -3.86201936  -0.07517988   5.13653145   5.13653145   7.91810064   9.68495616
+   3   0.34641016  -4.63079629   1.43419304   5.35284854   5.35284854   8.15540136   9.81498032
+   4   0.51961524  -5.25812925   3.25094996   5.69541211   5.69541211   8.51864422   9.65342396
+   5   0.69282032  -5.65198054   5.13706981   6.08466055   6.08466055   8.86784679   9.13706533
+   6   0.86602540  -5.78587376   6.28876728   6.28876728   6.28876728   8.83643550   8.83643550
diff --git a/tests/integrate/107_PW_OB_outputbands/refBANDS_1.dat b/tests/integrate/107_PW_OB_outputbands/refBANDS_1.dat
index ced64a27bf..0a1088ceb5 100644
--- a/tests/integrate/107_PW_OB_outputbands/refBANDS_1.dat
+++ b/tests/integrate/107_PW_OB_outputbands/refBANDS_1.dat
@@ -1,6 +1,6 @@
-1 0  -3.3870489 -0.79801666 5.0648903 5.0649025 7.841159 9.6059536
-2 0.17320508  -3.861988 -0.075180055 5.1365331 5.136532 7.9181141 9.6849637
-3 0.34641016  -4.6307948 1.43421 5.3529077 5.3528526 8.1553806 9.8149771
-4 0.51961524  -5.2581202 3.2509597 5.6954124 5.6954235 8.5188338 9.6533658
-5 0.69282032  -5.6519785 5.1370872 6.0846758 6.0846914 8.8679617 9.1370593
-6 0.8660254  -5.7858674 6.288817 6.2887755 6.2887815 8.8365307 8.8365071
+   1   0.00000000  -3.38704889  -0.79801666   5.06489026   5.06490253   7.84115900   9.60595365
+   2   0.17320508  -3.86198803  -0.07518005   5.13653307   5.13653200   7.91811409   9.68496369
+   3   0.34641016  -4.63079483   1.43421001   5.35290770   5.35285262   8.15538060   9.81497711
+   4   0.51961524  -5.25812025   3.25095973   5.69541240   5.69542353   8.51883375   9.65336577
+   5   0.69282032  -5.65197852   5.13708720   6.08467575   6.08469139   8.86796173   9.13705931
+   6   0.86602540  -5.78586742   6.28881704   6.28877545   6.28878152   8.83653074   8.83650706
diff --git a/tests/integrate/207_NO_KP_OB/refBANDS_1.dat b/tests/integrate/207_NO_KP_OB/refBANDS_1.dat
index c3cd3a0b6d..87a35be1b0 100644
--- a/tests/integrate/207_NO_KP_OB/refBANDS_1.dat
+++ b/tests/integrate/207_NO_KP_OB/refBANDS_1.dat
@@ -1,6 +1,6 @@
-1 0  -3.2007432 -0.55268317 5.3411505 5.3411505 8.2951616 10.292492
-2 0.17320508  -3.6638913 0.16302859 5.4122469 5.4122469 8.3831312 10.350824
-3 0.34641016  -4.4233335 1.6715127 5.6238068 5.6238068 8.6434268 10.399663
-4 0.51961524  -5.0448087 3.4994417 5.9592422 5.9592422 9.0296871 10.132257
-5 0.69282032  -5.4330707 5.3995429 6.3414385 6.3414385 9.3569963 9.5807977
-6 0.8660254  -5.5643351 6.5398122 6.5398122 6.5398122 9.280228 9.280228
+   1   0.00000000  -3.20074324  -0.55268317   5.34115051   5.34115051   8.29516160  10.29249200
+   2   0.17320508  -3.66389133   0.16302859   5.41224692   5.41224692   8.38313120  10.35082356
+   3   0.34641016  -4.42333349   1.67151268   5.62380682   5.62380682   8.64342675  10.39966330
+   4   0.51961524  -5.04480873   3.49944170   5.95924219   5.95924219   9.02968708  10.13225743
+   5   0.69282032  -5.43307067   5.39954292   6.34143853   6.34143853   9.35699627   9.58079775
+   6   0.86602540  -5.56433513   6.53981221   6.53981221   6.53981221   9.28022796   9.28022796

From 63875d6fee63b872aa6d737fa14c431cb5682157 Mon Sep 17 00:00:00 2001
From: Peng Xingliang <91927439+pxlxingliang@users.noreply.github.com>
Date: Fri, 19 Jan 2024 19:34:07 +0800
Subject: [PATCH 03/17] fix: bug in Autotest.sh when result.ref has no
 totaltimeref (#3523)

---
 tests/integrate/Autotest.sh | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/tests/integrate/Autotest.sh b/tests/integrate/Autotest.sh
index 5310998097..37908d1769 100755
--- a/tests/integrate/Autotest.sh
+++ b/tests/integrate/Autotest.sh
@@ -74,6 +74,11 @@ check_out(){
     # check every 'key' word
     #------------------------------------------------------
     for key in $properties; do
+    
+        if [ $key == "totaltimeref" ]; then
+            # echo "time=$cal ref=$ref"
+            break
+        fi
 
         #--------------------------------------------------
         # calculated value
@@ -91,11 +96,6 @@ check_out(){
         #--------------------------------------------------
         deviation=`awk 'BEGIN {x='$ref';y='$cal';printf "%.'$ca'f\n",x-y}'`
 
-        if [ $key == "totaltimeref" ]; then
-            # echo "time=$cal ref=$ref"
-            break
-        fi
-
 
         #--------------------------------------------------
         # If deviation < threshold, then the test passes,

From 4cd8d25f201e7221ce4c3bdbbdc8f50e22f5f09a Mon Sep 17 00:00:00 2001
From: Wenfei Li <38569667+wenfei-li@users.noreply.github.com>
Date: Fri, 19 Jan 2024 19:35:40 +0800
Subject: [PATCH 04/17] Fix : unit test of module_xc (#3524)

---
 deps/libpaw_interface                         |  2 +-
 .../module_xc/test/test_xc5.cpp               | 65 -------------------
 2 files changed, 1 insertion(+), 66 deletions(-)

diff --git a/deps/libpaw_interface b/deps/libpaw_interface
index 893cfe5b88..c211c0ab33 160000
--- a/deps/libpaw_interface
+++ b/deps/libpaw_interface
@@ -1 +1 @@
-Subproject commit 893cfe5b88c4b640b88a82335474d9f67d4c4cf6
+Subproject commit c211c0ab330adf3cc374f50ab3edee46b174e64c
diff --git a/source/module_hamilt_general/module_xc/test/test_xc5.cpp b/source/module_hamilt_general/module_xc/test/test_xc5.cpp
index bd2f87da8a..d9dfed1b20 100644
--- a/source/module_hamilt_general/module_xc/test/test_xc5.cpp
+++ b/source/module_hamilt_general/module_xc/test/test_xc5.cpp
@@ -25,9 +25,6 @@ class XCTest_VXC : public testing::Test
         double et2 = 0, vt2 = 0;
         ModuleBase::matrix v2;
 
-        double et4 = 0, vt4 = 0;
-        ModuleBase::matrix v4;
-
         void SetUp()
         {
             ModulePW::PW_Basis rhopw;
@@ -87,13 +84,6 @@ class XCTest_VXC : public testing::Test
             vt2 = std::get<1>(etxc_vtxc_v);
             v2  = std::get<2>(etxc_vtxc_v);
 
-            GlobalV::NSPIN = 4;
-            GlobalV::DOMAG = true;
-            etxc_vtxc_v
-                = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell);
-            et4 = std::get<0>(etxc_vtxc_v);
-            vt4 = std::get<1>(etxc_vtxc_v);
-            v4  = std::get<2>(etxc_vtxc_v);
         }
 };
 
@@ -121,29 +111,6 @@ TEST_F(XCTest_VXC, set_xc_type)
     EXPECT_NEAR(v2(1,3),-1.97506482,1.0e-8);
     EXPECT_NEAR(v2(1,4),-2.160374198,1.0e-8);
 
-    EXPECT_NEAR(et4,-27.40098253,1.0e-8);
-    EXPECT_NEAR(vt4,-35.81948838,1.0e-8);
-    EXPECT_NEAR(v4(0,0),0,1.0e-8);
-    EXPECT_NEAR(v4(0,1),-1.559604078,1.0e-8);
-    EXPECT_NEAR(v4(0,2),-1.920028447,1.0e-8);
-    EXPECT_NEAR(v4(0,3),-2.168396069,1.0e-8);
-    EXPECT_NEAR(v4(0,4),-2.36419592,1.0e-8);
-    EXPECT_NEAR(v4(1,0),0,1.0e-8);
-    EXPECT_NEAR(v4(1,1),-0.09308179605,1.0e-8);
-    EXPECT_NEAR(v4(1,2),-0.123132664,1.0e-8);
-    EXPECT_NEAR(v4(1,3),-0.144332804,1.0e-8);
-    EXPECT_NEAR(v4(1,4),-0.16127282,1.0e-8);
-    EXPECT_NEAR(v4(2,0),0,1.0e-8);
-    EXPECT_NEAR(v4(2,1),-0.9308179605,1.0e-8);
-    EXPECT_NEAR(v4(2,2),-1.23132664,1.0e-8);
-    EXPECT_NEAR(v4(2,3),-1.44332804,1.0e-8);
-    EXPECT_NEAR(v4(2,4),-1.6127282,1.0e-8);
-    EXPECT_NEAR(v4(3,0),0,1.0e-8);
-    EXPECT_NEAR(v4(3,1),-0.09308179605,1.0e-8);
-    EXPECT_NEAR(v4(3,2),-0.123132664,1.0e-8);
-    EXPECT_NEAR(v4(3,3),-0.144332804,1.0e-8);
-    EXPECT_NEAR(v4(3,4),-0.16127282,1.0e-8);
-
 }
 
 class XCTest_VXC_Libxc : public testing::Test
@@ -156,9 +123,6 @@ class XCTest_VXC_Libxc : public testing::Test
         double et2 = 0, vt2 = 0;
         ModuleBase::matrix v2;
 
-        double et4 = 0, vt4 = 0;
-        ModuleBase::matrix v4;
-
         void SetUp()
         {
             ModulePW::PW_Basis rhopw;
@@ -218,13 +182,6 @@ class XCTest_VXC_Libxc : public testing::Test
             vt2 = std::get<1>(etxc_vtxc_v);
             v2  = std::get<2>(etxc_vtxc_v);
 
-            GlobalV::NSPIN = 4;
-            GlobalV::DOMAG = true;
-            etxc_vtxc_v
-                = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell);
-            et4 = std::get<0>(etxc_vtxc_v);
-            vt4 = std::get<1>(etxc_vtxc_v);
-            v4  = std::get<2>(etxc_vtxc_v);
         }
 };
 
@@ -252,28 +209,6 @@ TEST_F(XCTest_VXC_Libxc, set_xc_type)
     EXPECT_NEAR(v2(1,3),-1.975058937,1.0e-8);
     EXPECT_NEAR(v2(1,4),-2.160368003,1.0e-8);
 
-    EXPECT_NEAR(et4,-27.28201062,1.0e-8);
-    EXPECT_NEAR(vt4,-35.98253991,1.0e-8);
-    EXPECT_NEAR(v4(0,0),0,1.0e-8);
-    EXPECT_NEAR(v4(0,1),-1.268278149,1.0e-8);
-    EXPECT_NEAR(v4(0,2),-1.598108222,1.0e-8);
-    EXPECT_NEAR(v4(0,3),-1.828079634,1.0e-8);
-    EXPECT_NEAR(v4(0,4),-2.010634115,1.0e-8);
-    EXPECT_NEAR(v4(1,0),0,1.0e-8);
-    EXPECT_NEAR(v4(1,1),-0.1255782493,1.0e-8);
-    EXPECT_NEAR(v4(1,2),-0.1582362929,1.0e-8);
-    EXPECT_NEAR(v4(1,3),-0.1810068558,1.0e-8);
-    EXPECT_NEAR(v4(1,4),-0.1990824429,1.0e-8);
-    EXPECT_NEAR(v4(2,0),0,1.0e-8);
-    EXPECT_NEAR(v4(2,1),-1.255782493,1.0e-8);
-    EXPECT_NEAR(v4(2,2),-1.582362929,1.0e-8);
-    EXPECT_NEAR(v4(2,3),-1.810068558,1.0e-8);
-    EXPECT_NEAR(v4(2,4),-1.990824429,1.0e-8);
-    EXPECT_NEAR(v4(3,0),0,1.0e-8);
-    EXPECT_NEAR(v4(3,1),-0.1255782493,1.0e-8);
-    EXPECT_NEAR(v4(3,2),-0.1582362929,1.0e-8);
-    EXPECT_NEAR(v4(3,3),-0.1810068558,1.0e-8);
-    EXPECT_NEAR(v4(3,4),-0.1990824429,1.0e-8);
 }
 
 class XCTest_VXC_meta : public testing::Test

From 600e1a7f071a2e720ce65ded04b727beccaea536 Mon Sep 17 00:00:00 2001
From: Zhao Tianqi <hongriTianqi@users.noreply.github.com>
Date: Mon, 22 Jan 2024 10:18:13 +0800
Subject: [PATCH 05/17] Fix: omit small magnetic moments to avoid numerical
 instability (#3530)

* update deltalambda

* avoid numerical error in orbMulP

* add constrain on Mi

* change case reference value
---
 .../module_deltaspin/cal_mw_helper.cpp        |  14 +-
 .../module_deltaspin/lambda_loop.cpp          |   4 +-
 .../204_NO_KP_NC_deltaspin/mulliken.txt.ref   | 168 +++++++++---------
 .../204_NO_KP_NC_deltaspin/result.ref         |   6 +-
 4 files changed, 94 insertions(+), 98 deletions(-)

diff --git a/source/module_hamilt_lcao/module_deltaspin/cal_mw_helper.cpp b/source/module_hamilt_lcao/module_deltaspin/cal_mw_helper.cpp
index 0c0595a9ce..bd0ad4ce3a 100644
--- a/source/module_hamilt_lcao/module_deltaspin/cal_mw_helper.cpp
+++ b/source/module_hamilt_lcao/module_deltaspin/cal_mw_helper.cpp
@@ -22,7 +22,7 @@ std::vector<std::vector<std::vector<double>>> SpinConstrain<std::complex<double>
                 AorbMulP[is][iat].resize(nw_it, 0.0);
                 for (int iw = 0; iw < nw_it; iw++)
                 {
-                    AorbMulP[is][iat][iw] = orbMulP(is, num);
+                    AorbMulP[is][iat][iw] = std::abs(orbMulP(is, num))< 1e-10 ? 0.0 : orbMulP(is, num);
                     num++;
                 }
             }
@@ -92,16 +92,10 @@ void SpinConstrain<std::complex<double>, psi::DEVICE_CPU>::calculate_MW(
             }
             else if (this->nspin_ == 4)
             {
-                this->Mi_[iat].x = total_charge_soc[1];
-                this->Mi_[iat].y = total_charge_soc[2];
-                this->Mi_[iat].z = total_charge_soc[3];
+                this->Mi_[iat].x = (std::abs(total_charge_soc[1]) < this->sc_thr_)? 0.0 : total_charge_soc[1];
+                this->Mi_[iat].y = (std::abs(total_charge_soc[2]) < this->sc_thr_)? 0.0 : total_charge_soc[2];
+                this->Mi_[iat].z = (std::abs(total_charge_soc[3]) < this->sc_thr_)? 0.0 : total_charge_soc[3];
             }
-            if (std::abs(this->Mi_[iat].x) < 1e-12)
-                this->Mi_[iat].x = 0.0;
-            if (std::abs(this->Mi_[iat].y) < 1e-12)
-                this->Mi_[iat].y = 0.0;
-            if (std::abs(this->Mi_[iat].z) < 1e-12)
-                this->Mi_[iat].z = 0.0;
         }
     }
 }
diff --git a/source/module_hamilt_lcao/module_deltaspin/lambda_loop.cpp b/source/module_hamilt_lcao/module_deltaspin/lambda_loop.cpp
index db64d5490d..845db88062 100644
--- a/source/module_hamilt_lcao/module_deltaspin/lambda_loop.cpp
+++ b/source/module_hamilt_lcao/module_deltaspin/lambda_loop.cpp
@@ -45,6 +45,7 @@ void SpinConstrain<std::complex<double>, psi::DEVICE_CPU>::run_lambda_loop(int o
         }
         else
         {
+            where_fill_scalar_else_2d(this->constrain_, 0, zero, delta_lambda, delta_lambda);
             add_scalar_multiply_2d(initial_lambda, delta_lambda, one, this->lambda_);
             this->cal_mw_from_lambda(i_step);
             new_spin = this->Mi_;
@@ -87,6 +88,7 @@ void SpinConstrain<std::complex<double>, psi::DEVICE_CPU>::run_lambda_loop(int o
         add_scalar_multiply_2d(dnu, search, alpha_trial, dnu);
         delta_lambda = dnu;
 
+        where_fill_scalar_else_2d(this->constrain_, 0, zero, delta_lambda, delta_lambda);
         add_scalar_multiply_2d(initial_lambda, delta_lambda, one, this->lambda_);
         this->cal_mw_from_lambda(i_step);
 
@@ -115,4 +117,4 @@ void SpinConstrain<std::complex<double>, psi::DEVICE_CPU>::run_lambda_loop(int o
         }
         alpha_trial = alpha_trial * pow(g, 0.7);
     }
-}
\ No newline at end of file
+}
diff --git a/tests/integrate/204_NO_KP_NC_deltaspin/mulliken.txt.ref b/tests/integrate/204_NO_KP_NC_deltaspin/mulliken.txt.ref
index dee68c8f40..bffad6b08a 100644
--- a/tests/integrate/204_NO_KP_NC_deltaspin/mulliken.txt.ref
+++ b/tests/integrate/204_NO_KP_NC_deltaspin/mulliken.txt.ref
@@ -3,92 +3,92 @@ CALCULATE THE MULLIkEN ANALYSIS FOR EACH ATOM
  Total charge:	32
 Decomposed Mulliken populations
 0                 Zeta of Fe                        Spin 1                        Spin 2                        Spin 3                        Spin 4
-s                        0                           1.317                       0.05552                       -0.2843                       0.02903
-  sum over m                                         1.317                       0.05552                       -0.2843                       0.02903
-s                        1                           1.726                      -0.01923                       0.09498                      0.005159
-  sum over m                                         1.726                      -0.01923                       0.09498                      0.005159
-s                        2                         0.03246                      -0.04333                        0.2148                      0.008137
-  sum over m                                       0.03246                      -0.04333                        0.2148                      0.008137
-s                        3                        -0.02921                      0.005194                      -0.02641                      0.001867
-  sum over m                                      -0.02921                      0.005194                      -0.02641                      0.001867
-  sum over m+zeta                                    3.046                     -0.001842                    -0.0009368                       0.04419
-pz                        0                           2.034                     -0.001185                      0.005932                     1.545e-06
-px                        0                           2.033                     -0.001283                      0.006419                     1.538e-06
-py                        0                           2.033                     -0.001188                      0.005944                     1.543e-06
-  sum over m                                           6.1                     -0.003656                       0.01829                     4.626e-06
-pz                        1                        -0.02622                     0.0005602                     -0.002791                             0
-px                        1                        -0.02639                     0.0006145                     -0.003054                             0
-py                        1                        -0.02603                     0.0005563                      -0.00277                             0
-  sum over m                                      -0.07864                      0.001731                     -0.008615                             0
-  sum over m+zeta                                    6.021                     -0.001925                       0.00968                     5.611e-06
-dz^2                        0                           1.964                     0.0008273                     -0.004131                     4.077e-06
-dxz                        0                           1.044                        0.1755                       -0.7507                      0.002258
-dyz                        0                          0.9544                        0.1768                       -0.7532                      0.002329
-dx^2-y^2                        0                           1.967                     0.0007523                     -0.003756                     3.978e-06
-dxy                        0                           1.055                        0.1751                       -0.7495                      0.002251
-  sum over m                                         6.984                         0.529                        -2.261                      0.006846
-dz^2                        1                         0.03863                    -0.0008699                      0.004363                     5.197e-06
-dxz                        1                        -0.03759                     -0.005346                       0.01936                    -0.0001322
-dyz                        1                        -0.03407                     -0.005734                       0.02118                    -0.0001342
-dx^2-y^2                        1                         0.03943                    -0.0009093                      0.004564                     5.691e-06
-dxy                        1                        -0.03787                     -0.005246                        0.0189                    -0.0001314
-  sum over m                                      -0.03146                      -0.01811                       0.06836                     -0.000387
-  sum over m+zeta                                    6.952                        0.5109                        -2.193                      0.006459
-fz^3                        0                       -0.007049                     0.0007578                     -0.003775                             0
-fxz^2                        0                       -0.002045                     0.0002638                     -0.001312                             0
-fyz^2                        0                       -0.002729                     0.0002912                     -0.001448                             0
-fzx^2-zy^2                        0                       6.273e-05                             0                     6.642e-06                             0
-fxyz                        0                       1.153e-05                     1.446e-06                    -5.675e-06                             0
-fx^3-3*xy^2                        0                        -0.00338                       0.00044                     -0.002189                             0
-f3yx^2-y^3                        0                        -0.00407                     0.0004646                     -0.002311                             0
-  sum over m                                       -0.0192                      0.002219                      -0.01103                     2.581e-06
-  sum over m+zeta                                  -0.0192                      0.002219                      -0.01103                     2.581e-06
-Total Charge on atom:  Fe                  16
-Total Magnetism on atom:  Fe                   (0.5093, -2.195, 0.05066)
+s                        0                           1.317                       0.06196                       -0.2625                      -0.07949
+  sum over m                                         1.317                       0.06196                       -0.2625                      -0.07949
+s                        1                           1.726                      -0.01809                       0.09886                      -0.01413
+  sum over m                                         1.726                      -0.01809                       0.09886                      -0.01413
+s                        2                         0.03246                      -0.04153                        0.2209                      -0.02228
+  sum over m                                       0.03246                      -0.04153                        0.2209                      -0.02228
+s                        3                        -0.02921                      0.005609                        -0.025                     -0.005114
+  sum over m                                      -0.02921                      0.005609                        -0.025                     -0.005114
+  sum over m+zeta                                    3.046                      0.007945                        0.0323                        -0.121
+pz                        0                           2.034                     -0.001186                      0.005932                    -3.981e-06
+px                        0                           2.033                     -0.001283                      0.006419                    -3.989e-06
+py                        0                           2.033                     -0.001188                      0.005944                    -3.979e-06
+  sum over m                                           6.1                     -0.003658                        0.0183                    -1.195e-05
+pz                        1                        -0.02621                     0.0005578                     -0.002789                             0
+px                        1                        -0.02639                     0.0006107                     -0.003054                             0
+py                        1                        -0.02603                     0.0005536                     -0.002768                             0
+  sum over m                                      -0.07863                      0.001722                     -0.008611                             0
+  sum over m+zeta                                    6.021                     -0.001936                      0.009684                    -1.277e-05
+dz^2                        0                           1.964                     0.0008269                     -0.004128                    -1.088e-05
+dxz                        0                           1.044                         0.156                       -0.7849                    -0.0003055
+dyz                        0                          0.9592                        0.1564                       -0.7869                    -0.0003096
+dx^2-y^2                        0                           1.967                      0.000752                     -0.003754                    -1.059e-05
+dxy                        0                           1.055                        0.1558                       -0.7835                    -0.0003047
+  sum over m                                         6.988                        0.4698                        -2.363                    -0.0009413
+dz^2                        1                         0.03863                    -0.0008716                      0.004365                    -1.357e-05
+dxz                        1                        -0.03708                     -0.004148                       0.02101                     1.956e-05
+dyz                        1                        -0.03373                     -0.004494                       0.02274                     1.968e-05
+dx^2-y^2                        1                         0.03943                    -0.0009117                      0.004566                    -1.471e-05
+dxy                        1                        -0.03733                     -0.004056                       0.02055                     1.945e-05
+  sum over m                                      -0.03008                      -0.01448                       0.07324                     3.041e-05
+  sum over m+zeta                                    6.958                        0.4553                         -2.29                    -0.0009109
+fz^3                        0                       -0.007044                     0.0007552                     -0.003776                    -1.406e-06
+fxz^2                        0                       -0.002046                     0.0002628                     -0.001314                             0
+fyz^2                        0                        -0.00273                       0.00029                      -0.00145                             0
+fzx^2-zy^2                        0                       5.811e-05                             0                     3.451e-06                             0
+fxyz                        0                        1.14e-05                     1.249e-06                    -6.306e-06                             0
+fx^3-3*xy^2                        0                       -0.003379                     0.0004381                      -0.00219                             0
+f3yx^2-y^3                        0                        -0.00407                     0.0004626                     -0.002313                             0
+  sum over m                                       -0.0192                      0.002209                      -0.01105                    -4.307e-06
+  sum over m+zeta                                  -0.0192                      0.002209                      -0.01105                    -4.307e-06
+Total Charge on atom:  Fe               16.01
+Total Magnetism on atom:  Fe                   (0.4635, -2.259, -0.1219)
 
 
 1                 Zeta of Fe                        Spin 1                        Spin 2                        Spin 3                        Spin 4
-s                        0                           1.275                       0.05341                       -0.2605                      -0.02903
-  sum over m                                         1.275                       0.05341                       -0.2605                      -0.02903
-s                        1                           1.755                      -0.01752                       0.08879                     -0.005156
-  sum over m                                         1.755                      -0.01752                       0.08879                     -0.005156
-s                        2                        -0.02898                       -0.0404                        0.2039                      -0.00813
-  sum over m                                      -0.02898                       -0.0404                        0.2039                      -0.00813
-s                        3                        -0.04711                      0.006367                      -0.03139                     -0.001874
-  sum over m                                      -0.04711                      0.006367                      -0.03139                     -0.001874
-  sum over m+zeta                                    2.954                      0.001862                     0.0008532                      -0.04419
-pz                        0                           2.032                     -0.001369                      0.006852                    -1.367e-06
-px                        0                           2.025                    -0.0009208                      0.004608                    -1.387e-06
-py                        0                           2.032                     -0.001332                      0.006666                    -1.366e-06
-  sum over m                                         6.089                     -0.003622                       0.01813                    -4.119e-06
-pz                        1                        -0.02528                     0.0005889                     -0.002889                             0
-px                        1                        -0.01606                     0.0001369                    -0.0006408                             0
-py                        1                        -0.02466                      0.000571                     -0.002802                             0
-  sum over m                                        -0.066                      0.001297                     -0.006331                     2.367e-06
-  sum over m+zeta                                    6.023                     -0.002325                       0.01179                    -1.753e-06
-dz^2                        0                           1.957                      0.001158                     -0.005774                    -3.913e-06
-dxz                        0                           1.097                        0.1724                       -0.7275                      0.002311
-dyz                        0                          0.9509                        0.1759                       -0.7475                      0.002269
-dx^2-y^2                        0                           1.947                      0.001654                     -0.008245                    -4.075e-06
-dxy                        0                           1.113                        0.1714                       -0.7227                      0.002304
-  sum over m                                         7.065                        0.5225                        -2.212                      0.006876
-dz^2                        1                         0.03925                     -0.001062                      0.005333                    -4.383e-06
-dxz                        1                         -0.0366                     -0.003947                       0.01263                    -0.0001213
-dyz                        1                        -0.03157                     -0.005197                       0.01856                    -0.0001267
-dx^2-y^2                        1                         0.04266                     -0.001394                      0.007002                    -4.206e-06
-dxy                        1                        -0.03743                     -0.003854                       0.01222                    -0.0001203
-  sum over m                                      -0.02369                      -0.01545                       0.05575                    -0.0003768
-  sum over m+zeta                                    7.041                        0.5071                        -2.156                      0.006499
-fz^3                        0                       -0.006614                     0.0007261                     -0.003596                             0
-fxz^2                        0                       -0.001954                     0.0002565                     -0.001276                             0
-fyz^2                        0                       -0.002684                     0.0002742                     -0.001366                             0
-fzx^2-zy^2                        0                        9.09e-05                      1.99e-05                    -8.018e-05                             0
-fxyz                        0                       2.062e-05                     4.102e-06                    -1.816e-05                             0
-fx^3-3*xy^2                        0                       -0.003203                     0.0004291                      -0.00213                             0
-f3yx^2-y^3                        0                       -0.003698                     0.0004635                     -0.002271                             0
-  sum over m                                      -0.01804                      0.002174                      -0.01074                             0
-  sum over m+zeta                                 -0.01804                      0.002174                      -0.01074                             0
-Total Charge on atom:  Fe                  16
-Total Magnetism on atom:  Fe                   (0.5088, -2.154, -0.03769)
+s                        0                           1.275                       0.04699                       -0.2823                       0.07949
+  sum over m                                         1.275                       0.04699                       -0.2823                       0.07949
+s                        1                           1.755                      -0.01866                       0.08491                       0.01412
+  sum over m                                         1.755                      -0.01866                       0.08491                       0.01412
+s                        2                        -0.02899                      -0.04221                        0.1978                       0.02226
+  sum over m                                      -0.02899                      -0.04221                        0.1978                       0.02226
+s                        3                        -0.04712                       0.00595                      -0.03281                      0.005133
+  sum over m                                      -0.04712                       0.00595                      -0.03281                      0.005133
+  sum over m+zeta                                    2.954                     -0.007928                      -0.03239                         0.121
+pz                        0                           2.032                     -0.001371                       0.00685                     3.967e-06
+px                        0                           2.025                    -0.0009218                      0.004606                     3.958e-06
+py                        0                           2.032                     -0.001333                      0.006664                     3.965e-06
+  sum over m                                         6.089                     -0.003626                       0.01812                     1.189e-05
+pz                        1                        -0.02529                     0.0005803                     -0.002904                             0
+px                        1                        -0.01606                     0.0001295                    -0.0006492                             0
+py                        1                        -0.02466                     0.0005625                     -0.002815                             0
+  sum over m                                      -0.06602                      0.001272                     -0.006367                             0
+  sum over m+zeta                                    6.023                     -0.002353                       0.01175                      1.25e-05
+dz^2                        0                           1.957                      0.001154                     -0.005778                     1.149e-05
+dxz                        0                           1.091                        0.1517                       -0.7637                    -8.462e-05
+dyz                        0                          0.9556                        0.1553                       -0.7815                    -8.443e-05
+dx^2-y^2                        0                           1.947                      0.001648                     -0.008249                     1.233e-05
+dxy                        0                           1.106                        0.1508                       -0.7591                    -8.432e-05
+  sum over m                                         7.056                        0.4606                        -2.318                    -0.0002295
+dz^2                        1                         0.03925                     -0.001067                      0.005328                     1.289e-05
+dxz                        1                        -0.03558                     -0.002824                       0.01439                       2.5e-06
+dyz                        1                        -0.03117                     -0.003962                        0.0201                     2.798e-06
+dx^2-y^2                        1                         0.04266                     -0.001401                      0.006997                      1.29e-05
+dxy                        1                        -0.03637                     -0.002747                       0.01401                     2.475e-06
+  sum over m                                      -0.02122                        -0.012                       0.06082                     3.356e-05
+  sum over m+zeta                                    7.035                        0.4486                        -2.257                     -0.000196
+fz^3                        0                       -0.006615                     0.0007206                     -0.003605                     1.352e-06
+fxz^2                        0                       -0.001955                     0.0002554                     -0.001278                             0
+fyz^2                        0                       -0.002684                     0.0002735                     -0.001368                             0
+fzx^2-zy^2                        0                       9.383e-05                      1.68e-05                    -8.473e-05                             0
+fxyz                        0                       2.053e-05                      3.66e-06                    -1.839e-05                             0
+fx^3-3*xy^2                        0                       -0.003204                     0.0004266                     -0.002134                             0
+f3yx^2-y^3                        0                       -0.003695                     0.0004558                     -0.002281                             0
+  sum over m                                      -0.01804                      0.002152                      -0.01077                     4.022e-06
+  sum over m+zeta                                 -0.01804                      0.002152                      -0.01077                     4.022e-06
+Total Charge on atom:  Fe               15.99
+Total Magnetism on atom:  Fe                   (0.4405, -2.289, 0.1208)
 
 
diff --git a/tests/integrate/204_NO_KP_NC_deltaspin/result.ref b/tests/integrate/204_NO_KP_NC_deltaspin/result.ref
index a5ca4b1941..8a17a1fada 100644
--- a/tests/integrate/204_NO_KP_NC_deltaspin/result.ref
+++ b/tests/integrate/204_NO_KP_NC_deltaspin/result.ref
@@ -1,4 +1,4 @@
-etotref -6844.685232778258
-etotperatomref -3422.3426163891
+etotref -6844.326716364628
+etotperatomref -3422.1633581823
 Compare_mulliken_pass 0
-totaltimeref 17.11
+totaltimeref 36.59

From 433d0efd1954b34d0aed98c39cf92fdb000f6b8c Mon Sep 17 00:00:00 2001
From: Chun Cai <amoycaic@gmail.com>
Date: Mon, 22 Jan 2024 10:19:29 +0800
Subject: [PATCH 06/17] Fix: fix multiple compiler warnings (#3515)

* Fix: add noreturn attribute to warning_quit

* Add type conversion

* fix string literal

* fix small number trunctuation

* Fix system call returned value not checked

* fix missing braket

* Refactor parameter_pool.cpp and parameter_pool.h

* remove duplicated return statements

* Change WARNING_QUIT occurances in tests

* Add warning message to help debug UT
---
 .../module_base/test/complexmatrix_test.cpp   |  22 +-
 .../module_base/test/inverse_matrix_test.cpp  |   2 +-
 source/module_base/test/math_ylmreal_test.cpp | 364 +++++++++---------
 source/module_base/tool_quit.h                |   8 +-
 .../test/sltk_atom_input_test.cpp             |   2 +-
 .../module_xc/test/test_xc.cpp                |   4 +-
 .../module_xc/test/test_xc1.cpp               |   2 +-
 .../module_xc/test/test_xc2.cpp               |  10 +-
 .../module_xc/test/test_xc4.cpp               |   4 +-
 .../module_xc/test/xc3_mock.h                 |  16 +-
 .../module_tddft/test/tddft_test.cpp          |   3 +-
 source/module_io/parameter_pool.cpp           |  33 +-
 source/module_io/parameter_pool.h             |   8 +-
 source/module_ri/LRI_CV_Tools.hpp             |  11 +-
 source/module_ri/exx_lip.cpp                  |  14 +-
 15 files changed, 248 insertions(+), 255 deletions(-)

diff --git a/source/module_base/test/complexmatrix_test.cpp b/source/module_base/test/complexmatrix_test.cpp
index 026aeb40de..0adc52363a 100644
--- a/source/module_base/test/complexmatrix_test.cpp
+++ b/source/module_base/test/complexmatrix_test.cpp
@@ -23,8 +23,8 @@
  *  - set_as_identity_matrix()
  *  - print():Output the elements of this complex matrix greater than threshold.
  *  - checkreal()
- * 
- * Tested relative functions 
+ *
+ * Tested relative functions
  *  - operator "+" "-" "*" between two ComplexMatrix
  *  - operator "*" between a ComplexMatrix and double or complex, and reverse.
  *  - trace()
@@ -35,13 +35,13 @@
  *  - conj()
  *  - scale_accumulate():
  *  - scaled_sum():
- * 
+ *
  */
 
 //a mock function of WARNING_QUIT, to avoid the uncorrected call by matrix.cpp at line 37.
 namespace ModuleBase
 {
-    void WARNING_QUIT(const std::string &file,const std::string &description) {return ;}
+    void WARNING_QUIT(const std::string &file,const std::string &description) {exit(1);}
 }
 
 inline void EXPECT_COMPLEX_EQ(const std::complex<double>& a,const std::complex<double>& b)
@@ -104,8 +104,8 @@ TEST_F(ComplexMatrixTest,ConstructorCM)
 
 TEST_F(ComplexMatrixTest,ConstructorCMrvalue)
 {
-    ModuleBase::ComplexMatrix cm2(cm22);        
-    ModuleBase::ComplexMatrix cm1(std::move(cm22)); 
+    ModuleBase::ComplexMatrix cm2(cm22);
+    ModuleBase::ComplexMatrix cm1(std::move(cm22));
     EXPECT_EQ(cm1.nr,cm2.nr);
     EXPECT_EQ(cm1.nc,cm2.nc);
     EXPECT_EQ(cm1.size,cm2.size);
@@ -338,15 +338,15 @@ TEST_F(ComplexMatrixTest,OperatorMultMatrix)
     EXPECT_EQ(cm33.nr,3);
     EXPECT_EQ(cm33.nc,3);
     EXPECT_EQ(cm33.size,9);
-    EXPECT_COMPLEX_EQ(cm33(0,0),std::complex<double>{-46.0,72.0  }); 
+    EXPECT_COMPLEX_EQ(cm33(0,0),std::complex<double>{-46.0,72.0  });
     EXPECT_COMPLEX_EQ(cm33(0,1),std::complex<double>{-46.0,118.0 });
     EXPECT_COMPLEX_EQ(cm33(0,2),std::complex<double>{-46.0,164.0 });
     EXPECT_COMPLEX_EQ(cm33(1,0),std::complex<double>{-54.0,84.0  });
     EXPECT_COMPLEX_EQ(cm33(1,1),std::complex<double>{-54.0,138.0 });
-    EXPECT_COMPLEX_EQ(cm33(1,2),std::complex<double>{-54.0,192.0 }); 
+    EXPECT_COMPLEX_EQ(cm33(1,2),std::complex<double>{-54.0,192.0 });
     EXPECT_COMPLEX_EQ(cm33(2,0),std::complex<double>{-62.0,96.0  });
     EXPECT_COMPLEX_EQ(cm33(2,1),std::complex<double>{-62.0,158.0 });
-    EXPECT_COMPLEX_EQ(cm33(2,2),std::complex<double>{-62.0,220.0 }); 
+    EXPECT_COMPLEX_EQ(cm33(2,2),std::complex<double>{-62.0,220.0 });
 
     EXPECT_DEATH(cm22 * cm32,"");
 }
@@ -525,7 +525,7 @@ TEST_F(ComplexMatrixTest,ScaleSumArray)
     cmout = new ModuleBase::ComplexMatrix*[2];
     cmin1 = new ModuleBase::ComplexMatrix*[2];
     cmin2 = new ModuleBase::ComplexMatrix*[2];
-   
+
     cmin1[0] = &cm1;
     cmin1[1] = &cm2;
     cmin2[0] = &cm3;
@@ -563,7 +563,7 @@ TEST_F(ComplexMatrixTest,print)
    EXPECT_THAT(output,testing::HasSubstr("(3,4)\t(4,5)\t"));
    ifs.close();
    remove("printtest1.log");
-// The condition of  std::abs(data)>threshold_abs && std::imag(data)) <= threshold_imag 
+// The condition of  std::abs(data)>threshold_abs && std::imag(data)) <= threshold_imag
    ofs.open("printtest2.log");
    cm22.print(ofs,1e-10,2);
    ofs.close();
diff --git a/source/module_base/test/inverse_matrix_test.cpp b/source/module_base/test/inverse_matrix_test.cpp
index df68f58a56..a871f906cd 100644
--- a/source/module_base/test/inverse_matrix_test.cpp
+++ b/source/module_base/test/inverse_matrix_test.cpp
@@ -19,7 +19,7 @@
 //a mock function of WARNING_QUIT, to avoid the uncorrected call by matrix.cpp at line 37.
 namespace ModuleBase
 {
-    void WARNING_QUIT(const std::string &file,const std::string &description) {return ;}
+    void WARNING_QUIT(const std::string &file,const std::string &description) {exit(1);}
 }
 
 TEST(InverseMatrixComplexTest, InverseMatrixComplex)
diff --git a/source/module_base/test/math_ylmreal_test.cpp b/source/module_base/test/math_ylmreal_test.cpp
index d5e7a504ed..13d0bd2b69 100644
--- a/source/module_base/test/math_ylmreal_test.cpp
+++ b/source/module_base/test/math_ylmreal_test.cpp
@@ -13,16 +13,16 @@
 ***********************************************/
 
 /**
- * For lmax <5 cases, the reference values are calculated by the formula from 
+ * For lmax <5 cases, the reference values are calculated by the formula from
  * https://formulasearchengine.com/wiki/Table_of_spherical_harmonics. Note, these
- * formula lack of the Condon–Shortley phase (-1)^m, and in this unit test, item 
+ * formula lack of the Condon–Shortley phase (-1)^m, and in this unit test, item
  * (-1)^m is multiplied.
  * For lmax >=5, the reference values are calculated by YlmReal::Ylm_Real.
  *
  * - Tested functions of class YlmReal
  *      - Ylm_Real
  *      - Ylm_Real2
- *      - rlylm 
+ *      - rlylm
  *      - YlmRealTemplate (double and float)
  *
  * - Tested functions of class Ylm
@@ -30,9 +30,9 @@
  *      - sph_harm
  *      - rl_sph_harm
  *      - grad_rl_sph_harm
- *      - equality_value_test: test the eqaulity of Ylm function between rl_sph_harm (spherical input) and  get_ylm_real (Cartesian input) 
+ *      - equality_value_test: test the eqaulity of Ylm function between rl_sph_harm (spherical input) and  get_ylm_real (Cartesian input)
  *      - equality_gradient_test:test the eqaulity of Ylm gradient function between grad_rl_sph_harm(spherical input) and  rlylm (Cartesian input)
- * 
+ *
  */
 
 
@@ -40,7 +40,7 @@
 //mock functions of WARNING_QUIT and WARNING
 namespace ModuleBase
 {
-    void WARNING_QUIT(const std::string &file,const std::string &description) {return ;}
+    void WARNING_QUIT(const std::string &file,const std::string &description) {exit(1);}
     void WARNING(const std::string &file,const std::string &description) {return ;}
 }
 
@@ -58,7 +58,7 @@ class YlmRealTest : public testing::Test
     ModuleBase::Vector3<double> *g; //vectors of the 4 points
     double *ref;        //reference of Ylm
     double *rly;        //Ylm
-    double (*rlgy)[3];  //the gradient of Ylm  
+    double (*rlgy)[3];  //the gradient of Ylm
     std::vector<double> rlyvector; //Ylm
     std::vector<std::vector<double>> rlgyvector; //the gradient of Ylm
 
@@ -91,101 +91,101 @@ class YlmRealTest : public testing::Test
     double y4m4(const double &x, const double &y, const double &z) {double r=norm(x,y,z); return 3./4.*sqrt(35./M_PI) * x*y*(x*x - y*y) / (r*r*r*r);}
 
     //the reference values are calculated by ModuleBase::Ylm::grad_rl_sph_harm
-    //1st dimension: example, 2nd dimension: Ylm, 3rd dimension: dx/dy/dz 
+    //1st dimension: example, 2nd dimension: Ylm, 3rd dimension: dx/dy/dz
     double rlgyref[4][64][3] = {
-        {   { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  4.88603e-01}, {-4.88603e-01,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00, -4.88603e-01,  0.00000e+00}, {-6.30783e-01,  0.00000e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00, -1.09255e+00}, 
-            { 0.00000e+00, -0.00000e+00,  0.00000e+00}, { 1.09255e+00,  0.00000e+00,  0.00000e+00}, {-0.00000e+00,  1.09255e+00, -0.00000e+00}, 
-            {-0.00000e+00,  0.00000e+00, -1.11953e+00}, { 1.37114e+00,  0.00000e+00, -0.00000e+00}, { 0.00000e+00,  4.57046e-01,  0.00000e+00}, 
-            { 0.00000e+00,  0.00000e+00,  1.44531e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, {-1.77013e+00,  0.00000e+00, -0.00000e+00}, 
-            { 0.00000e+00, -1.77013e+00,  0.00000e+00}, { 1.26943e+00,  0.00000e+00, -0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.00714e+00}, 
-            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, {-1.89235e+00,  0.00000e+00,  0.00000e+00}, {-0.00000e+00, -9.46175e-01,  0.00000e+00}, 
-            {-0.00000e+00,  0.00000e+00, -1.77013e+00}, { 0.00000e+00, -0.00000e+00,  0.00000e+00}, { 2.50334e+00,  0.00000e+00,  0.00000e+00}, 
-            {-0.00000e+00,  2.50334e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  1.75425e+00}, {-2.26473e+00,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00, -4.52947e-01,  0.00000e+00}, {-0.00000e+00,  0.00000e+00, -2.39677e+00}, {-0.00000e+00, -0.00000e+00,  0.00000e+00}, 
-            { 2.44619e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  1.46771e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.07566e+00}, 
-            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, {-3.28191e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -3.28191e+00,  0.00000e+00}, 
-            {-1.90708e+00,  0.00000e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00, -2.91311e+00}, { 0.00000e+00, -0.00000e+00,  0.00000e+00}, 
-            { 2.76362e+00,  0.00000e+00, -0.00000e+00}, {-0.00000e+00,  9.21205e-01,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.76362e+00}, 
-            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, {-3.02739e+00,  0.00000e+00,  0.00000e+00}, {-0.00000e+00, -2.01826e+00,  0.00000e+00}, 
-            {-0.00000e+00,  0.00000e+00, -2.36662e+00}, { 0.00000e+00, -0.00000e+00,  0.00000e+00}, { 4.09910e+00,  0.00000e+00,  0.00000e+00}, 
-            {-0.00000e+00,  4.09910e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00, -2.38995e+00}, { 3.16161e+00,  0.00000e+00, -0.00000e+00}, 
-            { 0.00000e+00,  4.51658e-01,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  3.31900e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            {-3.28564e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -1.40813e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00, -3.11349e+00}, 
-            {-0.00000e+00, -0.00000e+00,  0.00000e+00}, { 3.63241e+00,  0.00000e+00, -0.00000e+00}, { 0.00000e+00,  2.59458e+00,  0.00000e+00}, 
-            { 0.00000e+00,  0.00000e+00,  2.64596e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, {-4.95014e+00,  0.00000e+00, -0.00000e+00}, 
+        {   { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  4.88603e-01}, {-4.88603e-01,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00, -4.88603e-01,  0.00000e+00}, {-6.30783e-01,  0.00000e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00, -1.09255e+00},
+            { 0.00000e+00, -0.00000e+00,  0.00000e+00}, { 1.09255e+00,  0.00000e+00,  0.00000e+00}, {-0.00000e+00,  1.09255e+00, -0.00000e+00},
+            {-0.00000e+00,  0.00000e+00, -1.11953e+00}, { 1.37114e+00,  0.00000e+00, -0.00000e+00}, { 0.00000e+00,  4.57046e-01,  0.00000e+00},
+            { 0.00000e+00,  0.00000e+00,  1.44531e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, {-1.77013e+00,  0.00000e+00, -0.00000e+00},
+            { 0.00000e+00, -1.77013e+00,  0.00000e+00}, { 1.26943e+00,  0.00000e+00, -0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.00714e+00},
+            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, {-1.89235e+00,  0.00000e+00,  0.00000e+00}, {-0.00000e+00, -9.46175e-01,  0.00000e+00},
+            {-0.00000e+00,  0.00000e+00, -1.77013e+00}, { 0.00000e+00, -0.00000e+00,  0.00000e+00}, { 2.50334e+00,  0.00000e+00,  0.00000e+00},
+            {-0.00000e+00,  2.50334e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  1.75425e+00}, {-2.26473e+00,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00, -4.52947e-01,  0.00000e+00}, {-0.00000e+00,  0.00000e+00, -2.39677e+00}, {-0.00000e+00, -0.00000e+00,  0.00000e+00},
+            { 2.44619e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  1.46771e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.07566e+00},
+            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, {-3.28191e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -3.28191e+00,  0.00000e+00},
+            {-1.90708e+00,  0.00000e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00, -2.91311e+00}, { 0.00000e+00, -0.00000e+00,  0.00000e+00},
+            { 2.76362e+00,  0.00000e+00, -0.00000e+00}, {-0.00000e+00,  9.21205e-01,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.76362e+00},
+            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, {-3.02739e+00,  0.00000e+00,  0.00000e+00}, {-0.00000e+00, -2.01826e+00,  0.00000e+00},
+            {-0.00000e+00,  0.00000e+00, -2.36662e+00}, { 0.00000e+00, -0.00000e+00,  0.00000e+00}, { 4.09910e+00,  0.00000e+00,  0.00000e+00},
+            {-0.00000e+00,  4.09910e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00, -2.38995e+00}, { 3.16161e+00,  0.00000e+00, -0.00000e+00},
+            { 0.00000e+00,  4.51658e-01,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  3.31900e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00},
+            {-3.28564e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -1.40813e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00, -3.11349e+00},
+            {-0.00000e+00, -0.00000e+00,  0.00000e+00}, { 3.63241e+00,  0.00000e+00, -0.00000e+00}, { 0.00000e+00,  2.59458e+00,  0.00000e+00},
+            { 0.00000e+00,  0.00000e+00,  2.64596e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, {-4.95014e+00,  0.00000e+00, -0.00000e+00},
             { 0.00000e+00, -4.95014e+00,  0.00000e+00}
         },
         {
-            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  4.88603e-01}, {-4.88603e-01,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00, -4.88603e-01,  0.00000e+00}, { 0.00000e+00, -6.30783e-01,  0.00000e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00, -0.00000e+00, -1.09255e+00}, { 0.00000e+00, -1.09255e+00,  0.00000e+00}, { 1.09255e+00,  0.00000e+00, -0.00000e+00}, 
-            { 0.00000e+00, -0.00000e+00, -1.11953e+00}, { 4.57046e-01,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  1.37114e+00, -0.00000e+00}, 
-            { 0.00000e+00, -0.00000e+00, -1.44531e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 1.77013e+00,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00,  1.77013e+00,  0.00000e+00}, { 0.00000e+00,  1.26943e+00, -0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00,  0.00000e+00,  2.00714e+00}, { 0.00000e+00,  1.89235e+00, -0.00000e+00}, {-9.46175e-01,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  1.77013e+00}, { 0.00000e+00,  2.50334e+00, -0.00000e+00}, 
-            {-2.50334e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  1.75425e+00}, {-4.52947e-01,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00, -2.26473e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.39677e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            {-1.46771e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -2.44619e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.07566e+00}, 
-            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, {-3.28191e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -3.28191e+00,  0.00000e+00}, 
-            { 0.00000e+00, -1.90708e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -0.00000e+00, -2.91311e+00}, 
-            { 0.00000e+00, -2.76362e+00,  0.00000e+00}, { 9.21205e-01,  0.00000e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00, -0.00000e+00, -2.76362e+00}, { 0.00000e+00, -3.02739e+00,  0.00000e+00}, { 2.01826e+00,  0.00000e+00,  0.00000e+00}, 
-            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -0.00000e+00, -2.36662e+00}, { 0.00000e+00, -4.09910e+00,  0.00000e+00}, 
-            { 4.09910e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -0.00000e+00, -2.38995e+00}, { 4.51658e-01,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00,  3.16161e+00, -0.00000e+00}, { 0.00000e+00, -0.00000e+00, -3.31900e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            { 1.40813e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  3.28564e+00, -0.00000e+00}, { 0.00000e+00, -0.00000e+00, -3.11349e+00}, 
-            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 2.59458e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  3.63241e+00, -0.00000e+00}, 
-            { 0.00000e+00,  0.00000e+00, -2.64596e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 4.95014e+00,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00,  4.95014e+00, -0.00000e+00}           
+            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  4.88603e-01}, {-4.88603e-01,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00, -4.88603e-01,  0.00000e+00}, { 0.00000e+00, -6.30783e-01,  0.00000e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00, -0.00000e+00, -1.09255e+00}, { 0.00000e+00, -1.09255e+00,  0.00000e+00}, { 1.09255e+00,  0.00000e+00, -0.00000e+00},
+            { 0.00000e+00, -0.00000e+00, -1.11953e+00}, { 4.57046e-01,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  1.37114e+00, -0.00000e+00},
+            { 0.00000e+00, -0.00000e+00, -1.44531e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 1.77013e+00,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00,  1.77013e+00,  0.00000e+00}, { 0.00000e+00,  1.26943e+00, -0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00,  0.00000e+00,  2.00714e+00}, { 0.00000e+00,  1.89235e+00, -0.00000e+00}, {-9.46175e-01,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  1.77013e+00}, { 0.00000e+00,  2.50334e+00, -0.00000e+00},
+            {-2.50334e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  1.75425e+00}, {-4.52947e-01,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00, -2.26473e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.39677e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00},
+            {-1.46771e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -2.44619e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.07566e+00},
+            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, {-3.28191e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -3.28191e+00,  0.00000e+00},
+            { 0.00000e+00, -1.90708e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -0.00000e+00, -2.91311e+00},
+            { 0.00000e+00, -2.76362e+00,  0.00000e+00}, { 9.21205e-01,  0.00000e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00, -0.00000e+00, -2.76362e+00}, { 0.00000e+00, -3.02739e+00,  0.00000e+00}, { 2.01826e+00,  0.00000e+00,  0.00000e+00},
+            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -0.00000e+00, -2.36662e+00}, { 0.00000e+00, -4.09910e+00,  0.00000e+00},
+            { 4.09910e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -0.00000e+00, -2.38995e+00}, { 4.51658e-01,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00,  3.16161e+00, -0.00000e+00}, { 0.00000e+00, -0.00000e+00, -3.31900e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00},
+            { 1.40813e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  3.28564e+00, -0.00000e+00}, { 0.00000e+00, -0.00000e+00, -3.11349e+00},
+            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 2.59458e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  3.63241e+00, -0.00000e+00},
+            { 0.00000e+00,  0.00000e+00, -2.64596e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 4.95014e+00,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00,  4.95014e+00, -0.00000e+00}
         },
         {
-            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  4.88603e-01}, {-4.88603e-01,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00, -4.88603e-01,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  1.26157e+00}, {-1.09255e+00,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00, -1.09255e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.22045e-16}, {-0.00000e+00,  0.00000e+00, -0.00000e+00}, 
-            { 0.00000e+00,  0.00000e+00,  2.23906e+00}, {-1.82818e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -1.82818e+00,  0.00000e+00}, 
-            { 0.00000e+00,  0.00000e+00,  8.81212e-16}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, {-1.84324e-16,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00,  5.55112e-17,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  3.38514e+00}, {-2.67619e+00,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00, -2.67619e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.30756e-15}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            {-5.52973e-16,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  1.66533e-16,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  4.67801e+00}, {-3.62357e+00,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00, -3.62357e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  4.87108e-15}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            {-1.22267e-15,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  3.68219e-16,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, { 4.93038e-32,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -6.16298e-33,  0.00000e+00}, 
-            { 0.00000e+00,  0.00000e+00,  6.10264e+00}, {-4.66097e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -4.66097e+00,  0.00000e+00}, 
-            { 0.00000e+00,  0.00000e+00,  8.98664e-15}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, {-2.30221e-15,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00,  6.93334e-16,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            { 1.77767e-31,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -2.22209e-32,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  7.64784e+00}, {-5.78122e+00,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00, -5.78122e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  1.51096e-14}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            {-3.91011e-15,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  1.17757e-15,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, 
-            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, { 4.67737e-31,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -5.84671e-32,  0.00000e+00}, 
-            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 1.13319e-47,  0.00000e+00,  0.00000e+00}, 
+            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  4.88603e-01}, {-4.88603e-01,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00, -4.88603e-01,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  1.26157e+00}, {-1.09255e+00,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00, -1.09255e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.22045e-16}, {-0.00000e+00,  0.00000e+00, -0.00000e+00},
+            { 0.00000e+00,  0.00000e+00,  2.23906e+00}, {-1.82818e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -1.82818e+00,  0.00000e+00},
+            { 0.00000e+00,  0.00000e+00,  8.81212e-16}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, {-1.84324e-16,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00,  5.55112e-17,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  3.38514e+00}, {-2.67619e+00,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00, -2.67619e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  2.30756e-15}, {-0.00000e+00,  0.00000e+00,  0.00000e+00},
+            {-5.52973e-16,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  1.66533e-16,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00},
+            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  4.67801e+00}, {-3.62357e+00,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00, -3.62357e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  4.87108e-15}, {-0.00000e+00,  0.00000e+00,  0.00000e+00},
+            {-1.22267e-15,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  3.68219e-16,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00},
+            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, { 4.93038e-32,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -6.16298e-33,  0.00000e+00},
+            { 0.00000e+00,  0.00000e+00,  6.10264e+00}, {-4.66097e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -4.66097e+00,  0.00000e+00},
+            { 0.00000e+00,  0.00000e+00,  8.98664e-15}, {-0.00000e+00,  0.00000e+00,  0.00000e+00}, {-2.30221e-15,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00,  6.93334e-16,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, {-0.00000e+00,  0.00000e+00,  0.00000e+00},
+            { 1.77767e-31,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -2.22209e-32,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00},
+            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  7.64784e+00}, {-5.78122e+00,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00, -5.78122e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  1.51096e-14}, {-0.00000e+00,  0.00000e+00,  0.00000e+00},
+            {-3.91011e-15,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  1.17757e-15,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00},
+            {-0.00000e+00,  0.00000e+00,  0.00000e+00}, { 4.67737e-31,  0.00000e+00,  0.00000e+00}, { 0.00000e+00, -5.84671e-32,  0.00000e+00},
+            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 1.13319e-47,  0.00000e+00,  0.00000e+00},
             { 0.00000e+00, -1.41649e-48,  0.00000e+00}
         },
         {
-            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  4.88603e-01}, {-4.88603e-01,  0.00000e+00,  0.00000e+00}, 
-            { 0.00000e+00, -4.88603e-01,  0.00000e+00}, { 3.64183e-01,  3.64183e-01, -7.28366e-01}, { 6.30783e-01, -0.00000e+00,  6.30783e-01}, 
-            {-0.00000e+00,  6.30783e-01,  6.30783e-01}, {-6.30783e-01,  6.30783e-01, -1.66533e-16}, {-6.30783e-01, -6.30783e-01,  0.00000e+00}, 
-            {-7.46353e-01, -7.46353e-01,  0.00000e+00}, { 0.00000e+00,  3.04697e-01, -1.21879e+00}, { 3.04697e-01,  0.00000e+00, -1.21879e+00}, 
-            { 9.63537e-01, -9.63537e-01,  4.01253e-16}, { 9.63537e-01,  9.63537e-01,  9.63537e-01}, {-4.44089e-16,  1.18009e+00, -2.22045e-16}, 
-            {-1.18009e+00, -1.11022e-16,  0.00000e+00}, { 4.88603e-01,  4.88603e-01,  1.30294e+00}, {-1.03006e+00, -7.72548e-01,  7.72548e-01}, 
-            {-7.72548e-01, -1.03006e+00,  7.72548e-01}, {-7.28366e-01,  7.28366e-01, -5.25363e-16}, {-3.64183e-01, -3.64183e-01, -2.18510e+00}, 
-            { 7.69185e-16, -2.04397e+00, -6.81324e-01}, { 2.04397e+00,  1.92296e-16,  6.81324e-01}, { 9.63537e-01,  9.63537e-01, -1.44756e-16}, 
-            {-9.63537e-01,  9.63537e-01, -5.55112e-17}, { 5.19779e-01,  5.19779e-01, -1.81923e+00}, { 1.40917e+00,  8.05238e-01,  8.05238e-01}, 
-            { 8.05238e-01,  1.40917e+00,  8.05238e-01}, { 0.00000e+00, -4.44089e-16,  3.24739e-16}, {-1.06523e+00, -1.06523e+00,  2.13046e+00}, 
-            {-2.17439e-01,  1.73951e+00,  1.73951e+00}, {-1.73951e+00,  2.17439e-01, -1.73951e+00}, {-1.84503e+00, -1.84503e+00, -9.22517e-01}, 
-            { 1.84503e+00, -1.84503e+00,  6.58625e-16}, { 1.45863e+00,  1.11022e-15,  0.00000e+00}, {-8.88178e-16,  1.45863e+00,  0.00000e+00}, 
-            {-1.46807e+00, -1.46807e+00,  5.87227e-01}, {-4.48502e-01, -3.36617e-16, -2.24251e+00}, {-3.36617e-16, -4.48502e-01, -2.24251e+00}, 
-            { 7.09144e-01, -7.09144e-01,  1.87222e-16}, { 2.12743e+00,  2.12743e+00, -9.38779e-16}, { 7.09144e-01, -5.11006e-16, -2.12743e+00}, 
-            { 1.02201e-15, -7.09144e-01,  2.12743e+00}, { 1.81260e+00,  1.81260e+00,  2.58943e+00}, {-2.07154e+00,  2.07154e+00, -1.66969e-15}, 
-            {-3.03637e+00, -2.31111e-15, -6.07275e-01}, { 1.84889e-15, -3.03637e+00, -6.07275e-01}, { 1.05183e+00, -1.05183e+00,  5.77778e-17}, 
-            { 1.05183e+00,  1.05183e+00,  4.03986e-17}, { 1.27464e+00,  1.27464e+00,  1.69952e+00}, {-1.28472e+00, -1.20442e+00,  1.92707e+00}, 
-            {-1.20442e+00, -1.28472e+00,  1.92707e+00}, {-8.52285e-01,  8.52285e-01, -6.74704e-16}, {-1.50789e+00, -1.50789e+00, -2.95022e+00}, 
-            {-1.11260e+00, -2.08612e+00,  9.27164e-01}, { 2.08612e+00,  1.11260e+00, -9.27164e-01}, {-3.07506e-01, -3.07506e-01, -3.69007e+00}, 
-            { 1.23002e+00, -1.23002e+00,  2.28018e-15}, { 3.69007e+00, -1.53753e-01,  1.84503e+00}, {-1.53753e-01,  3.69007e+00,  1.84503e+00}, 
-            {-2.35197e+00,  2.35197e+00, -8.00513e-16}, {-2.35197e+00, -2.35197e+00, -7.83988e-01}, { 1.37903e-15, -1.46671e+00,  9.77875e-17}, 
+            { 0.00000e+00,  0.00000e+00,  0.00000e+00}, { 0.00000e+00,  0.00000e+00,  4.88603e-01}, {-4.88603e-01,  0.00000e+00,  0.00000e+00},
+            { 0.00000e+00, -4.88603e-01,  0.00000e+00}, { 3.64183e-01,  3.64183e-01, -7.28366e-01}, { 6.30783e-01, -0.00000e+00,  6.30783e-01},
+            {-0.00000e+00,  6.30783e-01,  6.30783e-01}, {-6.30783e-01,  6.30783e-01, -1.66533e-16}, {-6.30783e-01, -6.30783e-01,  0.00000e+00},
+            {-7.46353e-01, -7.46353e-01,  0.00000e+00}, { 0.00000e+00,  3.04697e-01, -1.21879e+00}, { 3.04697e-01,  0.00000e+00, -1.21879e+00},
+            { 9.63537e-01, -9.63537e-01,  4.01253e-16}, { 9.63537e-01,  9.63537e-01,  9.63537e-01}, {-4.44089e-16,  1.18009e+00, -2.22045e-16},
+            {-1.18009e+00, -1.11022e-16,  0.00000e+00}, { 4.88603e-01,  4.88603e-01,  1.30294e+00}, {-1.03006e+00, -7.72548e-01,  7.72548e-01},
+            {-7.72548e-01, -1.03006e+00,  7.72548e-01}, {-7.28366e-01,  7.28366e-01, -5.25363e-16}, {-3.64183e-01, -3.64183e-01, -2.18510e+00},
+            { 7.69185e-16, -2.04397e+00, -6.81324e-01}, { 2.04397e+00,  1.92296e-16,  6.81324e-01}, { 9.63537e-01,  9.63537e-01, -1.44756e-16},
+            {-9.63537e-01,  9.63537e-01, -5.55112e-17}, { 5.19779e-01,  5.19779e-01, -1.81923e+00}, { 1.40917e+00,  8.05238e-01,  8.05238e-01},
+            { 8.05238e-01,  1.40917e+00,  8.05238e-01}, { 0.00000e+00, -4.44089e-16,  3.24739e-16}, {-1.06523e+00, -1.06523e+00,  2.13046e+00},
+            {-2.17439e-01,  1.73951e+00,  1.73951e+00}, {-1.73951e+00,  2.17439e-01, -1.73951e+00}, {-1.84503e+00, -1.84503e+00, -9.22517e-01},
+            { 1.84503e+00, -1.84503e+00,  6.58625e-16}, { 1.45863e+00,  1.11022e-15,  0.00000e+00}, {-8.88178e-16,  1.45863e+00,  0.00000e+00},
+            {-1.46807e+00, -1.46807e+00,  5.87227e-01}, {-4.48502e-01, -3.36617e-16, -2.24251e+00}, {-3.36617e-16, -4.48502e-01, -2.24251e+00},
+            { 7.09144e-01, -7.09144e-01,  1.87222e-16}, { 2.12743e+00,  2.12743e+00, -9.38779e-16}, { 7.09144e-01, -5.11006e-16, -2.12743e+00},
+            { 1.02201e-15, -7.09144e-01,  2.12743e+00}, { 1.81260e+00,  1.81260e+00,  2.58943e+00}, {-2.07154e+00,  2.07154e+00, -1.66969e-15},
+            {-3.03637e+00, -2.31111e-15, -6.07275e-01}, { 1.84889e-15, -3.03637e+00, -6.07275e-01}, { 1.05183e+00, -1.05183e+00,  5.77778e-17},
+            { 1.05183e+00,  1.05183e+00,  4.03986e-17}, { 1.27464e+00,  1.27464e+00,  1.69952e+00}, {-1.28472e+00, -1.20442e+00,  1.92707e+00},
+            {-1.20442e+00, -1.28472e+00,  1.92707e+00}, {-8.52285e-01,  8.52285e-01, -6.74704e-16}, {-1.50789e+00, -1.50789e+00, -2.95022e+00},
+            {-1.11260e+00, -2.08612e+00,  9.27164e-01}, { 2.08612e+00,  1.11260e+00, -9.27164e-01}, {-3.07506e-01, -3.07506e-01, -3.69007e+00},
+            { 1.23002e+00, -1.23002e+00,  2.28018e-15}, { 3.69007e+00, -1.53753e-01,  1.84503e+00}, {-1.53753e-01,  3.69007e+00,  1.84503e+00},
+            {-2.35197e+00,  2.35197e+00, -8.00513e-16}, {-2.35197e+00, -2.35197e+00, -7.83988e-01}, { 1.37903e-15, -1.46671e+00,  9.77875e-17},
             { 1.46671e+00,  1.14919e-15,  1.34475e-16}
         }
     };
@@ -206,71 +206,71 @@ class YlmRealTest : public testing::Test
         rlgy = new double[nylm][3];
         rlgyvector.resize(nylm,std::vector<double>(3));
         ref = new double[64*4]{
-            y00(g[0].x, g[0].y, g[0].z),  y00(g[1].x, g[1].y, g[1].z),  y00(g[2].x, g[2].y, g[2].z),  y00(g[3].x, g[3].y, g[3].z),  
-            y10(g[0].x, g[0].y, g[0].z),  y10(g[1].x, g[1].y, g[1].z),  y10(g[2].x, g[2].y, g[2].z),  y10(g[3].x, g[3].y, g[3].z),  
-            y11(g[0].x, g[0].y, g[0].z),  y11(g[1].x, g[1].y, g[1].z),  y11(g[2].x, g[2].y, g[2].z),  y11(g[3].x, g[3].y, g[3].z),  
-            y1m1(g[0].x, g[0].y, g[0].z), y1m1(g[1].x, g[1].y, g[1].z), y1m1(g[2].x, g[2].y, g[2].z), y1m1(g[3].x, g[3].y, g[3].z), 
-            y20(g[0].x, g[0].y, g[0].z),  y20(g[1].x, g[1].y, g[1].z),  y20(g[2].x, g[2].y, g[2].z),  y20(g[3].x, g[3].y, g[3].z),  
-            y21(g[0].x, g[0].y, g[0].z),  y21(g[1].x, g[1].y, g[1].z),  y21(g[2].x, g[2].y, g[2].z),  y21(g[3].x, g[3].y, g[3].z),  
-            y2m1(g[0].x, g[0].y, g[0].z), y2m1(g[1].x, g[1].y, g[1].z), y2m1(g[2].x, g[2].y, g[2].z), y2m1(g[3].x, g[3].y, g[3].z),                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               
-            y22(g[0].x, g[0].y, g[0].z),  y22(g[1].x, g[1].y, g[1].z),  y22(g[2].x, g[2].y, g[2].z),  y22(g[3].x, g[3].y, g[3].z),  
-            y2m2(g[0].x, g[0].y, g[0].z), y2m2(g[1].x, g[1].y, g[1].z), y2m2(g[2].x, g[2].y, g[2].z), y2m2(g[3].x, g[3].y, g[3].z), 
-            y30(g[0].x, g[0].y, g[0].z),  y30(g[1].x, g[1].y, g[1].z),  y30(g[2].x, g[2].y, g[2].z),  y30(g[3].x, g[3].y, g[3].z),  
-            y31(g[0].x, g[0].y, g[0].z),  y31(g[1].x, g[1].y, g[1].z),  y31(g[2].x, g[2].y, g[2].z),  y31(g[3].x, g[3].y, g[3].z),  
-            y3m1(g[0].x, g[0].y, g[0].z), y3m1(g[1].x, g[1].y, g[1].z), y3m1(g[2].x, g[2].y, g[2].z), y3m1(g[3].x, g[3].y, g[3].z), 
-            y32(g[0].x, g[0].y, g[0].z),  y32(g[1].x, g[1].y, g[1].z),  y32(g[2].x, g[2].y, g[2].z),  y32(g[3].x, g[3].y, g[3].z),  
-            y3m2(g[0].x, g[0].y, g[0].z), y3m2(g[1].x, g[1].y, g[1].z), y3m2(g[2].x, g[2].y, g[2].z), y3m2(g[3].x, g[3].y, g[3].z), 
-            y33(g[0].x, g[0].y, g[0].z),  y33(g[1].x, g[1].y, g[1].z),  y33(g[2].x, g[2].y, g[2].z),  y33(g[3].x, g[3].y, g[3].z),  
-            y3m3(g[0].x, g[0].y, g[0].z), y3m3(g[1].x, g[1].y, g[1].z), y3m3(g[2].x, g[2].y, g[2].z), y3m3(g[3].x, g[3].y, g[3].z), 
-            y40(g[0].x, g[0].y, g[0].z),  y40(g[1].x, g[1].y, g[1].z),  y40(g[2].x, g[2].y, g[2].z),  y40(g[3].x, g[3].y, g[3].z),  
-            y41(g[0].x, g[0].y, g[0].z),  y41(g[1].x, g[1].y, g[1].z),  y41(g[2].x, g[2].y, g[2].z),  y41(g[3].x, g[3].y, g[3].z),  
-            y4m1(g[0].x, g[0].y, g[0].z), y4m1(g[1].x, g[1].y, g[1].z), y4m1(g[2].x, g[2].y, g[2].z), y4m1(g[3].x, g[3].y, g[3].z), 
-            y42(g[0].x, g[0].y, g[0].z),  y42(g[1].x, g[1].y, g[1].z),  y42(g[2].x, g[2].y, g[2].z),  y42(g[3].x, g[3].y, g[3].z),  
-            y4m2(g[0].x, g[0].y, g[0].z), y4m2(g[1].x, g[1].y, g[1].z), y4m2(g[2].x, g[2].y, g[2].z), y4m2(g[3].x, g[3].y, g[3].z), 
-            y43(g[0].x, g[0].y, g[0].z),  y43(g[1].x, g[1].y, g[1].z),  y43(g[2].x, g[2].y, g[2].z),  y43(g[3].x, g[3].y, g[3].z),  
-            y4m3(g[0].x, g[0].y, g[0].z), y4m3(g[1].x, g[1].y, g[1].z), y4m3(g[2].x, g[2].y, g[2].z), y4m3(g[3].x, g[3].y, g[3].z), 
-            y44(g[0].x, g[0].y, g[0].z),  y44(g[1].x, g[1].y, g[1].z),  y44(g[2].x, g[2].y, g[2].z),  y44(g[3].x, g[3].y, g[3].z),  
-            y4m4(g[0].x, g[0].y, g[0].z), y4m4(g[1].x, g[1].y, g[1].z), y4m4(g[2].x, g[2].y, g[2].z), y4m4(g[3].x, g[3].y, g[3].z), 
-              0.000000000000000,    0.000000000000000,    0.935602579627389,    0.090028400200397, 
-             -0.452946651195697,   -0.000000000000000,   -0.000000000000000,   -0.348678494661834, 
-             -0.000000000000000,   -0.452946651195697,   -0.000000000000000,   -0.348678494661834, 
-             -0.000000000000000,    0.000000000000000,    0.000000000000000,   -0.000000000000000, 
-             -0.000000000000000,   -0.000000000000000,    0.000000000000000,   -0.000000000000000, 
-              0.489238299435250,    0.000000000000000,   -0.000000000000000,   -0.376615818502422, 
-              0.000000000000000,   -0.489238299435250,   -0.000000000000000,    0.376615818502422, 
-              0.000000000000000,    0.000000000000000,    0.000000000000000,    0.532615198330370, 
-              0.000000000000000,    0.000000000000000,    0.000000000000000,   -0.000000000000000, 
-             -0.656382056840170,   -0.000000000000000,   -0.000000000000000,   -0.168427714314628, 
-             -0.000000000000000,   -0.656382056840170,   -0.000000000000000,   -0.168427714314628, 
-             -0.317846011338142,   -0.317846011338142,    1.017107236282055,    0.226023830284901, 
-             -0.000000000000000,   -0.000000000000000,   -0.000000000000000,    0.258942827786103, 
-             -0.000000000000000,   -0.000000000000000,   -0.000000000000000,    0.258942827786103, 
-              0.460602629757462,   -0.460602629757462,    0.000000000000000,   -0.000000000000000, 
-              0.000000000000000,    0.000000000000000,    0.000000000000000,   -0.409424559784410, 
-             -0.000000000000000,   -0.000000000000000,   -0.000000000000000,    0.136474853261470, 
-             -0.000000000000000,    0.000000000000000,   -0.000000000000000,   -0.136474853261470, 
-             -0.504564900728724,   -0.504564900728724,    0.000000000000000,   -0.598002845308118, 
-             -0.000000000000000,   -0.000000000000000,    0.000000000000000,    0.000000000000000, 
-             -0.000000000000000,   -0.000000000000000,   -0.000000000000000,    0.350610246256556, 
-             -0.000000000000000,   -0.000000000000000,   -0.000000000000000,    0.350610246256556, 
-              0.683184105191914,   -0.683184105191914,    0.000000000000000,   -0.000000000000000, 
-              0.000000000000000,    0.000000000000000,    0.000000000000000,   -0.202424920056864, 
-              0.000000000000000,    0.000000000000000,    1.092548430592079,   -0.350435072502801, 
-              0.451658037912587,    0.000000000000000,   -0.000000000000000,    0.046358202625865, 
-              0.000000000000000,    0.451658037912587,   -0.000000000000000,    0.046358202625865, 
-              0.000000000000000,   -0.000000000000000,    0.000000000000000,    0.000000000000000, 
-              0.000000000000000,    0.000000000000000,    0.000000000000000,    0.492067081245654, 
-             -0.469376801586882,   -0.000000000000000,   -0.000000000000000,    0.187354445356332, 
-             -0.000000000000000,    0.469376801586882,   -0.000000000000000,   -0.187354445356332, 
-              0.000000000000000,    0.000000000000000,    0.000000000000000,    0.355076798886913, 
-              0.000000000000000,    0.000000000000000,    0.000000000000000,   -0.000000000000000, 
-              0.518915578720260,    0.000000000000000,   -0.000000000000000,   -0.443845998608641, 
-              0.000000000000000,    0.518915578720260,   -0.000000000000000,   -0.443845998608641, 
-              0.000000000000000,   -0.000000000000000,    0.000000000000000,    0.000000000000000, 
-              0.000000000000000,    0.000000000000000,    0.000000000000000,    0.452635881587108, 
-             -0.707162732524596,    0.000000000000000,   -0.000000000000000,    0.120972027847095, 
-             -0.000000000000000,    0.707162732524596,   -0.000000000000000,   -0.120972027847095  
-         } ; 
+            y00(g[0].x, g[0].y, g[0].z),  y00(g[1].x, g[1].y, g[1].z),  y00(g[2].x, g[2].y, g[2].z),  y00(g[3].x, g[3].y, g[3].z),
+            y10(g[0].x, g[0].y, g[0].z),  y10(g[1].x, g[1].y, g[1].z),  y10(g[2].x, g[2].y, g[2].z),  y10(g[3].x, g[3].y, g[3].z),
+            y11(g[0].x, g[0].y, g[0].z),  y11(g[1].x, g[1].y, g[1].z),  y11(g[2].x, g[2].y, g[2].z),  y11(g[3].x, g[3].y, g[3].z),
+            y1m1(g[0].x, g[0].y, g[0].z), y1m1(g[1].x, g[1].y, g[1].z), y1m1(g[2].x, g[2].y, g[2].z), y1m1(g[3].x, g[3].y, g[3].z),
+            y20(g[0].x, g[0].y, g[0].z),  y20(g[1].x, g[1].y, g[1].z),  y20(g[2].x, g[2].y, g[2].z),  y20(g[3].x, g[3].y, g[3].z),
+            y21(g[0].x, g[0].y, g[0].z),  y21(g[1].x, g[1].y, g[1].z),  y21(g[2].x, g[2].y, g[2].z),  y21(g[3].x, g[3].y, g[3].z),
+            y2m1(g[0].x, g[0].y, g[0].z), y2m1(g[1].x, g[1].y, g[1].z), y2m1(g[2].x, g[2].y, g[2].z), y2m1(g[3].x, g[3].y, g[3].z),
+            y22(g[0].x, g[0].y, g[0].z),  y22(g[1].x, g[1].y, g[1].z),  y22(g[2].x, g[2].y, g[2].z),  y22(g[3].x, g[3].y, g[3].z),
+            y2m2(g[0].x, g[0].y, g[0].z), y2m2(g[1].x, g[1].y, g[1].z), y2m2(g[2].x, g[2].y, g[2].z), y2m2(g[3].x, g[3].y, g[3].z),
+            y30(g[0].x, g[0].y, g[0].z),  y30(g[1].x, g[1].y, g[1].z),  y30(g[2].x, g[2].y, g[2].z),  y30(g[3].x, g[3].y, g[3].z),
+            y31(g[0].x, g[0].y, g[0].z),  y31(g[1].x, g[1].y, g[1].z),  y31(g[2].x, g[2].y, g[2].z),  y31(g[3].x, g[3].y, g[3].z),
+            y3m1(g[0].x, g[0].y, g[0].z), y3m1(g[1].x, g[1].y, g[1].z), y3m1(g[2].x, g[2].y, g[2].z), y3m1(g[3].x, g[3].y, g[3].z),
+            y32(g[0].x, g[0].y, g[0].z),  y32(g[1].x, g[1].y, g[1].z),  y32(g[2].x, g[2].y, g[2].z),  y32(g[3].x, g[3].y, g[3].z),
+            y3m2(g[0].x, g[0].y, g[0].z), y3m2(g[1].x, g[1].y, g[1].z), y3m2(g[2].x, g[2].y, g[2].z), y3m2(g[3].x, g[3].y, g[3].z),
+            y33(g[0].x, g[0].y, g[0].z),  y33(g[1].x, g[1].y, g[1].z),  y33(g[2].x, g[2].y, g[2].z),  y33(g[3].x, g[3].y, g[3].z),
+            y3m3(g[0].x, g[0].y, g[0].z), y3m3(g[1].x, g[1].y, g[1].z), y3m3(g[2].x, g[2].y, g[2].z), y3m3(g[3].x, g[3].y, g[3].z),
+            y40(g[0].x, g[0].y, g[0].z),  y40(g[1].x, g[1].y, g[1].z),  y40(g[2].x, g[2].y, g[2].z),  y40(g[3].x, g[3].y, g[3].z),
+            y41(g[0].x, g[0].y, g[0].z),  y41(g[1].x, g[1].y, g[1].z),  y41(g[2].x, g[2].y, g[2].z),  y41(g[3].x, g[3].y, g[3].z),
+            y4m1(g[0].x, g[0].y, g[0].z), y4m1(g[1].x, g[1].y, g[1].z), y4m1(g[2].x, g[2].y, g[2].z), y4m1(g[3].x, g[3].y, g[3].z),
+            y42(g[0].x, g[0].y, g[0].z),  y42(g[1].x, g[1].y, g[1].z),  y42(g[2].x, g[2].y, g[2].z),  y42(g[3].x, g[3].y, g[3].z),
+            y4m2(g[0].x, g[0].y, g[0].z), y4m2(g[1].x, g[1].y, g[1].z), y4m2(g[2].x, g[2].y, g[2].z), y4m2(g[3].x, g[3].y, g[3].z),
+            y43(g[0].x, g[0].y, g[0].z),  y43(g[1].x, g[1].y, g[1].z),  y43(g[2].x, g[2].y, g[2].z),  y43(g[3].x, g[3].y, g[3].z),
+            y4m3(g[0].x, g[0].y, g[0].z), y4m3(g[1].x, g[1].y, g[1].z), y4m3(g[2].x, g[2].y, g[2].z), y4m3(g[3].x, g[3].y, g[3].z),
+            y44(g[0].x, g[0].y, g[0].z),  y44(g[1].x, g[1].y, g[1].z),  y44(g[2].x, g[2].y, g[2].z),  y44(g[3].x, g[3].y, g[3].z),
+            y4m4(g[0].x, g[0].y, g[0].z), y4m4(g[1].x, g[1].y, g[1].z), y4m4(g[2].x, g[2].y, g[2].z), y4m4(g[3].x, g[3].y, g[3].z),
+              0.000000000000000,    0.000000000000000,    0.935602579627389,    0.090028400200397,
+             -0.452946651195697,   -0.000000000000000,   -0.000000000000000,   -0.348678494661834,
+             -0.000000000000000,   -0.452946651195697,   -0.000000000000000,   -0.348678494661834,
+             -0.000000000000000,    0.000000000000000,    0.000000000000000,   -0.000000000000000,
+             -0.000000000000000,   -0.000000000000000,    0.000000000000000,   -0.000000000000000,
+              0.489238299435250,    0.000000000000000,   -0.000000000000000,   -0.376615818502422,
+              0.000000000000000,   -0.489238299435250,   -0.000000000000000,    0.376615818502422,
+              0.000000000000000,    0.000000000000000,    0.000000000000000,    0.532615198330370,
+              0.000000000000000,    0.000000000000000,    0.000000000000000,   -0.000000000000000,
+             -0.656382056840170,   -0.000000000000000,   -0.000000000000000,   -0.168427714314628,
+             -0.000000000000000,   -0.656382056840170,   -0.000000000000000,   -0.168427714314628,
+             -0.317846011338142,   -0.317846011338142,    1.017107236282055,    0.226023830284901,
+             -0.000000000000000,   -0.000000000000000,   -0.000000000000000,    0.258942827786103,
+             -0.000000000000000,   -0.000000000000000,   -0.000000000000000,    0.258942827786103,
+              0.460602629757462,   -0.460602629757462,    0.000000000000000,   -0.000000000000000,
+              0.000000000000000,    0.000000000000000,    0.000000000000000,   -0.409424559784410,
+             -0.000000000000000,   -0.000000000000000,   -0.000000000000000,    0.136474853261470,
+             -0.000000000000000,    0.000000000000000,   -0.000000000000000,   -0.136474853261470,
+             -0.504564900728724,   -0.504564900728724,    0.000000000000000,   -0.598002845308118,
+             -0.000000000000000,   -0.000000000000000,    0.000000000000000,    0.000000000000000,
+             -0.000000000000000,   -0.000000000000000,   -0.000000000000000,    0.350610246256556,
+             -0.000000000000000,   -0.000000000000000,   -0.000000000000000,    0.350610246256556,
+              0.683184105191914,   -0.683184105191914,    0.000000000000000,   -0.000000000000000,
+              0.000000000000000,    0.000000000000000,    0.000000000000000,   -0.202424920056864,
+              0.000000000000000,    0.000000000000000,    1.092548430592079,   -0.350435072502801,
+              0.451658037912587,    0.000000000000000,   -0.000000000000000,    0.046358202625865,
+              0.000000000000000,    0.451658037912587,   -0.000000000000000,    0.046358202625865,
+              0.000000000000000,   -0.000000000000000,    0.000000000000000,    0.000000000000000,
+              0.000000000000000,    0.000000000000000,    0.000000000000000,    0.492067081245654,
+             -0.469376801586882,   -0.000000000000000,   -0.000000000000000,    0.187354445356332,
+             -0.000000000000000,    0.469376801586882,   -0.000000000000000,   -0.187354445356332,
+              0.000000000000000,    0.000000000000000,    0.000000000000000,    0.355076798886913,
+              0.000000000000000,    0.000000000000000,    0.000000000000000,   -0.000000000000000,
+              0.518915578720260,    0.000000000000000,   -0.000000000000000,   -0.443845998608641,
+              0.000000000000000,    0.518915578720260,   -0.000000000000000,   -0.443845998608641,
+              0.000000000000000,   -0.000000000000000,    0.000000000000000,    0.000000000000000,
+              0.000000000000000,    0.000000000000000,    0.000000000000000,    0.452635881587108,
+             -0.707162732524596,    0.000000000000000,   -0.000000000000000,    0.120972027847095,
+             -0.000000000000000,    0.707162732524596,   -0.000000000000000,   -0.120972027847095
+         } ;
     }
 
     void TearDown()
@@ -293,11 +293,11 @@ TEST_F(YlmRealTest,YlmReal)
     ModuleBase::YlmReal::Ylm_Real(nylm,ng,g,ylm);
     for(int i=0;i<nylm;++i)
     {
-        for(int j=0;j<ng;++j) 
+        for(int j=0;j<ng;++j)
         {
             EXPECT_NEAR(ylm(i,j),ref[i*ng+j],doublethreshold)  << "Ylm[" << i << "], example " << j << " not pass";
         }
-    } 
+    }
 }
 
 TEST_F(YlmRealTest,YlmRealTemplate)
@@ -318,7 +318,7 @@ TEST_F(YlmRealTest,gradYlmReal)
     ModuleBase::YlmReal::grad_Ylm_Real(nylm,ng,g,ylm,dylm[0],dylm[1],dylm[2]);
     for(int i=0;i<nylm;++i)
     {
-        for(int j=0;j<ng;++j) 
+        for(int j=0;j<ng;++j)
         {
             EXPECT_NEAR(ylm(i,j),ref[i*ng+j],doublethreshold)  << "Ylm[" << i << "], example " << j << " not pass";
         }
@@ -328,7 +328,7 @@ TEST_F(YlmRealTest,gradYlmReal)
     double step = 1e-7;
     for(int id = 0 ; id < 3 ; ++id)
     {
-        for(int j=0;j<ng;++j) 
+        for(int j=0;j<ng;++j)
         {
             ModuleBase::Vector3<double> gplus = g[j];
             ModuleBase::Vector3<double> gminus = g[j];
@@ -352,16 +352,16 @@ TEST_F(YlmRealTest,YlmReal2)
     ModuleBase::YlmReal::Ylm_Real2(nylm,ng,g,ylm);
     for(int i=0;i<nylm;++i)
     {
-        for(int j=0;j<ng;++j) 
+        for(int j=0;j<ng;++j)
         {
             EXPECT_NEAR(ylm(i,j),ref[i*ng+j],doublethreshold) << "Ylm[" << i << "], example " << j << " not pass";
         }
-    } 
+    }
 }
 
 
 TEST_F(YlmRealTest,YlmRealRlylm)
-{    
+{
     for(int j=0;j<ng;++j)
     {
         ModuleBase::YlmReal::rlylm(lmax,g[j].x,g[j].y,g[j].z,rly);
@@ -374,7 +374,7 @@ TEST_F(YlmRealTest,YlmRealRlylm)
 
 
 TEST_F(YlmRealTest,YlmGetYlmReal)
-{    
+{
     for(int j=0;j<ng;++j)
     {
         ModuleBase::Ylm::get_ylm_real(lmax+1,g[j],rly);
@@ -386,7 +386,7 @@ TEST_F(YlmRealTest,YlmGetYlmReal)
 }
 
 TEST_F(YlmRealTest,YlmSphHarm)
-{    
+{
     ModuleBase::Ylm::set_coefficients ();
     for(int j=0;j<ng;++j)
     {
@@ -395,13 +395,13 @@ TEST_F(YlmRealTest,YlmSphHarm)
         for(int i=0;i<nylm;++i)
         {
             EXPECT_NEAR(rlyvector[i],ref[i*ng+j],doublethreshold)  << "Ylm[" << i << "], example " << j << " not pass";
-            
+
         }
     }
 }
 
 TEST_F(YlmRealTest,YlmRlSphHarm)
-{    
+{
     ModuleBase::Ylm::set_coefficients ();
     for(int j=0;j<ng;++j)
     {
@@ -410,13 +410,13 @@ TEST_F(YlmRealTest,YlmRlSphHarm)
         for(int i=0;i<nylm;++i)
         {
             EXPECT_NEAR(rlyvector[i],ref[i*ng+j],doublethreshold)  << "Ylm[" << i << "], example " << j << " not pass";
-            
+
         }
     }
 }
 //used to be test1 in ylm.h
 TEST_F(YlmRealTest,YlmGradRlSphHarm)
-{    
+{
     ModuleBase::Ylm::set_coefficients ();
     for(int j=0;j<ng;++j)
     {
@@ -426,7 +426,7 @@ TEST_F(YlmRealTest,YlmGradRlSphHarm)
         {
             EXPECT_NEAR(rlyvector[i],ref[i*ng+j],doublethreshold)  << "Ylm[" << i << "], example " << j << " not pass";
             for(int k=0;k<3;++k) {EXPECT_NEAR(rlgyvector[i][k],rlgyref[j][i][k],1e-5);}
-            
+
         }
     }
 }
@@ -435,7 +435,7 @@ TEST_F(YlmRealTest,YlmGradRlSphHarm)
 TEST_F(YlmRealTest, equality_value_test)
 {
 
-    
+
     ModuleBase::Vector3<double> R (20.0, 0.0, 0.0);
 	const double xdr = R.x/R.norm();
 	const double ydr = R.y/R.norm();
@@ -444,17 +444,17 @@ TEST_F(YlmRealTest, equality_value_test)
 	const double rl = std::pow( R.norm(), L);
 	//std::cout << " rl=" << rl << std::endl;
 	ModuleBase::Ylm::set_coefficients();
-	
+
 	int nu = 100;
-	
+
 	// Peize Lin change rlya 2016-08-26
 	std::vector<double> rlya;
 	double rlyb[400];
 	ModuleBase::Ylm::ZEROS( rlyb, 400);
-	
+
 	ModuleBase::Ylm::rl_sph_harm(L, xdr, ydr, zdr, rlya);
 	ModuleBase::Ylm::get_ylm_real(L+1, R, rlyb);
-	
+
 	for (int i=0; i < nu; i++)
 	{
 		double diff = fabs(rlya[i]-rlyb[i]);
@@ -467,21 +467,21 @@ TEST_F(YlmRealTest, equality_value_test)
 TEST_F(YlmRealTest, equality_gradient_test)
 {
 
-    
+
     ModuleBase::Vector3<double> R (0.1,-0.2,0.5);
 	ModuleBase::Ylm::set_coefficients();
-	
+
 	//int nu = 100;
 
 	std::vector<double> rlya;
 	double rlyb[400];
-	
+
 	std::vector<std::vector<double>> grlya;
 	double grlyb[400][3];
-	
+
 	ModuleBase::Ylm::grad_rl_sph_harm (9, R.x, R.y, R.z, rlya, grlya);
 	ModuleBase::Ylm::rlylm (10, R.x, R.y, R.z, rlyb, grlyb);
-	
+
 	for (int i = 0; i < 100; i++)
 	{
 		double diffx = fabs(grlya[i][2]-grlyb[i][2]);
diff --git a/source/module_base/tool_quit.h b/source/module_base/tool_quit.h
index eafaf673cc..f944696d5a 100644
--- a/source/module_base/tool_quit.h
+++ b/source/module_base/tool_quit.h
@@ -33,13 +33,13 @@ void WARNING(const std::string &file, const std::string &description);
  * @brief Close .log files and exit
  *
  */
-void QUIT(void);
+[[noreturn]] void QUIT(void);
 
 /**
  * @brief Close .log files and exit
  *
  */
-void QUIT(int ret);
+[[noreturn]] void QUIT(int ret);
 
 /**
  * @brief Combine the functions of WARNING and QUIT
@@ -47,7 +47,7 @@ void QUIT(int ret);
  * @param file The file where warning happens
  * @param description The warning information
  */
-void WARNING_QUIT(const std::string &file, const std::string &description);
+[[noreturn]] void WARNING_QUIT(const std::string& file, const std::string& description);
 
 /**
  * @brief Combine the functions of WARNING and QUIT
@@ -55,7 +55,7 @@ void WARNING_QUIT(const std::string &file, const std::string &description);
  * @param file The file where warning happens
  * @param description The warning information
  */
-void WARNING_QUIT(const std::string &file, const std::string &description, int ret);
+[[noreturn]] void WARNING_QUIT(const std::string& file, const std::string& description, int ret);
 
 /**
  * @brief Check, if true, WARNING_QUIT
diff --git a/source/module_cell/module_neighbor/test/sltk_atom_input_test.cpp b/source/module_cell/module_neighbor/test/sltk_atom_input_test.cpp
index 617674256a..bb447bca4c 100644
--- a/source/module_cell/module_neighbor/test/sltk_atom_input_test.cpp
+++ b/source/module_cell/module_neighbor/test/sltk_atom_input_test.cpp
@@ -223,7 +223,7 @@ TEST_F(SltkAtomInputTest, ConstructorNoExpand)
     GlobalV::test_grid = 1;
     // this is a bug if radius is too small
     // because the expand_flag will be false!
-    radius = 1e-1000;
+    radius = 0;
     Atom_input Atom_inp(ofs, *ucell, ucell->nat, ucell->ntype, pbc, radius, test_atom_in);
     EXPECT_FALSE(Atom_inp.getExpandFlag());
     // call set_FAtom and Load_atom
diff --git a/source/module_hamilt_general/module_xc/test/test_xc.cpp b/source/module_hamilt_general/module_xc/test/test_xc.cpp
index 558556b66b..a770a88458 100644
--- a/source/module_hamilt_general/module_xc/test/test_xc.cpp
+++ b/source/module_hamilt_general/module_xc/test/test_xc.cpp
@@ -11,7 +11,7 @@
 
 namespace ModuleBase
 {
-    void WARNING_QUIT(const std::string &file,const std::string &description) {return ;}
+    void WARNING_QUIT(const std::string &file,const std::string &description) {exit(1);}
 }
 
 namespace GlobalV
@@ -95,7 +95,7 @@ class XCTest_PBEsol : public testing::Test
                 e_gga.push_back(e);
                 v1_gga.push_back(v1);
                 v2_gga.push_back(v2);
-            }                                           
+            }
         }
 };
 
diff --git a/source/module_hamilt_general/module_xc/test/test_xc1.cpp b/source/module_hamilt_general/module_xc/test/test_xc1.cpp
index 8e7a451e71..bc5c439630 100644
--- a/source/module_hamilt_general/module_xc/test/test_xc1.cpp
+++ b/source/module_hamilt_general/module_xc/test/test_xc1.cpp
@@ -12,7 +12,7 @@
 
 namespace ModuleBase
 {
-    void WARNING_QUIT(const std::string &file,const std::string &description) {return ;}
+    void WARNING_QUIT(const std::string &file,const std::string &description) {exit(1);}
 }
 
 namespace GlobalV
diff --git a/source/module_hamilt_general/module_xc/test/test_xc2.cpp b/source/module_hamilt_general/module_xc/test/test_xc2.cpp
index 4b1b7e888e..5bf75a3c68 100644
--- a/source/module_hamilt_general/module_xc/test/test_xc2.cpp
+++ b/source/module_hamilt_general/module_xc/test/test_xc2.cpp
@@ -11,7 +11,7 @@
 
 namespace ModuleBase
 {
-    void WARNING_QUIT(const std::string &file,const std::string &description) {return ;}
+    void WARNING_QUIT(const std::string &file,const std::string &description) {exit(1);}
 }
 
 namespace GlobalV
@@ -202,7 +202,7 @@ class XCTest_PZ_SPN : public testing::Test
                 e_lda.push_back(e);
                 v1_lda.push_back(v1);
                 v2_lda.push_back(v2);
-            }          
+            }
         }
 };
 
@@ -238,7 +238,7 @@ class XCTest_SLATER1_SPN : public testing::Test
                 e_lda.push_back(e);
                 v1_lda.push_back(v1);
                 v2_lda.push_back(v2);
-            }          
+            }
         }
 };
 
@@ -273,7 +273,7 @@ class XCTest_SLATER_RXC_SPN : public testing::Test
                 e_lda.push_back(e);
                 v1_lda.push_back(v1);
                 v2_lda.push_back(v2);
-            }           
+            }
         }
 };
 
@@ -310,7 +310,7 @@ class XCTest_P86_SPN : public testing::Test
                 v1_gga.push_back(v1);
                 v2_gga.push_back(v2);
                 v3_gga.push_back(v3);
-            }         
+            }
         }
 };
 
diff --git a/source/module_hamilt_general/module_xc/test/test_xc4.cpp b/source/module_hamilt_general/module_xc/test/test_xc4.cpp
index 114c817b0f..b4c8b70093 100644
--- a/source/module_hamilt_general/module_xc/test/test_xc4.cpp
+++ b/source/module_hamilt_general/module_xc/test/test_xc4.cpp
@@ -11,7 +11,7 @@
 
 namespace ModuleBase
 {
-    void WARNING_QUIT(const std::string &file,const std::string &description) {return ;}
+    void WARNING_QUIT(const std::string &file,const std::string &description) {exit(1);}
 }
 
 namespace GlobalV
@@ -49,7 +49,7 @@ class XCTest_SCAN : public testing::Test
                 v2_.push_back(v2);
                 v3_.push_back(v3);
             }
-        }                       
+        }
 };
 
 TEST_F(XCTest_SCAN, set_xc_type)
diff --git a/source/module_hamilt_general/module_xc/test/xc3_mock.h b/source/module_hamilt_general/module_xc/test/xc3_mock.h
index 628937adfe..da7f1e6f08 100644
--- a/source/module_hamilt_general/module_xc/test/xc3_mock.h
+++ b/source/module_hamilt_general/module_xc/test/xc3_mock.h
@@ -75,7 +75,7 @@ namespace ModulePW
         return x;
     }
 
-    
+
     template <typename FPTYPE, typename Device>
     void PW_Basis_K::real_to_recip(const Device* ctx,
                        const std::complex<FPTYPE>* in,
@@ -115,7 +115,7 @@ namespace ModulePW
                                                                      const int ik,
                                                                      const bool add,
                                                                      const double factor) const;
-#if __CUDA || __ROCM                                                                     
+#if __CUDA || __ROCM
     template void PW_Basis_K::real_to_recip<double, psi::DEVICE_GPU>(const psi::DEVICE_GPU* ctx,
                                                                      const std::complex<double>* in,
                                                                      std::complex<double>* out,
@@ -129,7 +129,7 @@ namespace ModulePW
                                                                      const int ik,
                                                                      const bool add,
                                                                      const double factor) const;
-#endif 
+#endif
 
     FFT::FFT(){};
     FFT::~FFT(){};
@@ -144,9 +144,13 @@ namespace ModulePW
 
 namespace ModuleBase
 {
-    void WARNING_QUIT(const std::string &file,const std::string &description) {return ;}
+    void WARNING_QUIT(const std::string &file,const std::string &description)
+    {
+        std::cout << " " << file <<"  warning : "<< description<<std::endl;
+        exit(1);
+    }
     void WARNING(const std::string &file,const std::string &description) {};
-    
+
     void Matrix3::Identity(){};
 
     IntArray::IntArray(int,int){};
@@ -228,4 +232,4 @@ namespace Parallel_Reduce
     template void reduce_pool<float>(float& object);
     template void reduce_pool<float>(float* object, const int n);
     template void reduce_pool<double>(double* object, const int n);
-}
\ No newline at end of file
+}
diff --git a/source/module_hamilt_lcao/module_tddft/test/tddft_test.cpp b/source/module_hamilt_lcao/module_tddft/test/tddft_test.cpp
index fedb46a976..a55ad59681 100644
--- a/source/module_hamilt_lcao/module_tddft/test/tddft_test.cpp
+++ b/source/module_hamilt_lcao/module_tddft/test/tddft_test.cpp
@@ -28,7 +28,8 @@ void MPIInit()
     npcol = 1;
     Cblacs_pinfo(&myrank, &mysize);
     Cblacs_get(-1, 0, &ictxt);
-    Cblacs_gridinit(&ictxt, "Row", nprow, npcol);
+    char order[] = "Row";
+    Cblacs_gridinit(&ictxt, order, nprow, npcol);
     Cblacs_gridinfo(ictxt, &nprow, &npcol, &myprow, &mypcol);
 }
 
diff --git a/source/module_io/parameter_pool.cpp b/source/module_io/parameter_pool.cpp
index 9a6fbed3f4..cf3e611ab3 100644
--- a/source/module_io/parameter_pool.cpp
+++ b/source/module_io/parameter_pool.cpp
@@ -69,7 +69,7 @@ int count_ntype(const std::string& fn)
  * @param input_value_path parameter default value file path
  * @param input_value_path parameter input value file path
  */
-bool Init(const std::string& default_type_path,
+void Init(const std::string& default_type_path,
           const std::string& default_value_path,
           const std::string& input_value_path)
 {
@@ -103,10 +103,8 @@ void strtolower(char* sa, char* sb)
  * @brief Reads the default parameters from the specified file and saves them to the global variable
  *        default_parametes_type
  * @param fn Specifies the path to the file
- * @return true Read successfully
- * @return false Read failure
  */
-bool default_parametes_reader(const std::string& fn, std::map<std::string, std::string>& default_parametes_type)
+void default_parametes_reader(const std::string& fn, std::map<std::string, std::string>& default_parametes_type)
 {
     std::ifstream inputFile(fn.c_str());
     if (inputFile.is_open())
@@ -122,28 +120,24 @@ bool default_parametes_reader(const std::string& fn, std::map<std::string, std::
     }
     else
     {
-        std::cout << "Cannot open file !" << std::endl;
+        ModuleBase::WARNING_QUIT("Input", "Cannot open file" + fn);
     }
 }
 /**
  * @brief This function is used to read the input parameter file and store it as a key-value pair
  * @param fn Enter the path to the parameter file
  */
-bool input_parameters_get(const std::string& fn, std::map<std::string, InputParameter>& input)
+void input_parameters_get(const std::string& fn, std::map<std::string, InputParameter>& input)
 {
-    // The module title information is displayed
     ModuleBase::TITLE("Input", "Read");
-    // If it is not the primary node, return false
     if (GlobalV::MY_RANK != 0)
-        return false;
+        return;
 
     // Open the input parameter file
     std::ifstream ifs(fn.c_str(), std::ios::in); // "in_datas/input_parameters"
-    // If the opening fails, an error message is printed and false is returned
     if (!ifs)
     {
-        std::cout << " Can't find the INPUT file." << std::endl;
-        return false;
+        ModuleBase::WARNING_QUIT("Input", "Can't find the INPUT file at " + fn);
     }
     ifs.clear();
     ifs.seekg(0);
@@ -166,8 +160,7 @@ bool input_parameters_get(const std::string& fn, std::map<std::string, InputPara
     // If ierr is 0, the word "INPUT_PARAMETERS" is not found, and an error message is printed with false
     if (ierr == 0)
     {
-        std::cout << " Error parameter list." << std::endl;
-        return false; // return error : false
+        ModuleBase::WARNING_QUIT("Input", "INPUT_PARAMETERS statement not found.");
     }
     ifs.rdstate();
 
@@ -274,15 +267,11 @@ bool input_parameters_get(const std::string& fn, std::map<std::string, InputPara
         }
         else if (ifs.bad() != 0)
         {
-            std::cout << " Bad input parameters. " << std::endl;
-            return false;
+            ModuleBase::WARNING_QUIT("Input", "Bad input parameters.");
         }
         else if (ifs.fail() != 0)
         {
-            std::cout << " word = " << word << std::endl;
-            std::cout << " Fail to read parameters. " << std::endl;
-            ifs.clear();
-            return false;
+            ModuleBase::WARNING_QUIT("Input", "Fail to read parameters: word = " + std::string(word));
         }
         else if (ifs.good() == 0)
         {
@@ -306,11 +295,9 @@ bool input_parameters_get(const std::string& fn, std::map<std::string, InputPara
     {
         ModuleBase::WARNING_QUIT("Input", "The ntype in INPUT is not equal to the ntype counted in STRU, check it.");
     }
-
-    return true;
 }
 
-bool input_parameters_set(std::map<std::string, InputParameter> input_parameters)
+void input_parameters_set(std::map<std::string, InputParameter> input_parameters)
 {
     if (input_parameters.count("nupdown") != 0)
     {
diff --git a/source/module_io/parameter_pool.h b/source/module_io/parameter_pool.h
index 83baedd036..bd4ae575dd 100644
--- a/source/module_io/parameter_pool.h
+++ b/source/module_io/parameter_pool.h
@@ -241,12 +241,12 @@ class InputParameter
         }
     }
 };
-bool Init(const std::string& default_type_path,
+void Init(const std::string& default_type_path,
           const std::string& default_value_path,
           const std::string& input_value_path);
-bool default_parametes_reader(const std::string& fn, std::map<std::string, std::string>& default_parametes_type);
-bool input_parameters_get(const std::string& fn, std::map<std::string, InputParameter>& input);
-bool input_parameters_set(std::map<std::string, InputParameter> input_parameters);
+void default_parametes_reader(const std::string& fn, std::map<std::string, std::string>& default_parametes_type);
+void input_parameters_get(const std::string& fn, std::map<std::string, InputParameter>& input);
+void input_parameters_set(std::map<std::string, InputParameter> input_parameters);
 
 extern std::map<std::string, InputParameter> input_parameters;
 extern std::map<std::string, std::string> default_parametes_type;
diff --git a/source/module_ri/LRI_CV_Tools.hpp b/source/module_ri/LRI_CV_Tools.hpp
index 532e7104fb..8ad95c3715 100644
--- a/source/module_ri/LRI_CV_Tools.hpp
+++ b/source/module_ri/LRI_CV_Tools.hpp
@@ -250,11 +250,10 @@ LRI_CV_Tools::cal_latvec_range(const double &rcut_times)
 	const ModuleBase::Vector3<double> proj = ModuleBase::Mathzone::latvec_projection(
 		std::array<ModuleBase::Vector3<double>,3>{GlobalC::ucell.a1, GlobalC::ucell.a2, GlobalC::ucell.a3});
 	const ModuleBase::Vector3<double> latvec_times = Rcut_max * rcut_times / (proj * GlobalC::ucell.lat0);
-	const ModuleBase::Vector3<Tcell> latvec_times_ceil = 
-		{std::ceil(latvec_times.x),
-		 std::ceil(latvec_times.y),
-		 std::ceil(latvec_times.z)};
-	const ModuleBase::Vector3<Tcell> period = 2 * latvec_times_ceil + ModuleBase::Vector3<Tcell>{1,1,1};
+    const ModuleBase::Vector3<Tcell> latvec_times_ceil = {static_cast<Tcell>(std::ceil(latvec_times.x)),
+                                                          static_cast<Tcell>(std::ceil(latvec_times.y)),
+                                                          static_cast<Tcell>(std::ceil(latvec_times.z))};
+    const ModuleBase::Vector3<Tcell> period = 2 * latvec_times_ceil + ModuleBase::Vector3<Tcell>{1,1,1};
 	return std::array<Tcell,3>{period.x, period.y, period.z};
 }
 
@@ -308,7 +307,7 @@ LRI_CV_Tools::get_dCVws(
 				const Abfs::Vector3_Order<double> R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*GlobalC::ucell.latvec);
 				dCVws[it0][it1][R_delta][ix] = dCVs_B.second;
 			}
-		}		
+		}
 	}
 	return dCVws;
 }
diff --git a/source/module_ri/exx_lip.cpp b/source/module_ri/exx_lip.cpp
index 2f685be5f1..0c4211d890 100644
--- a/source/module_ri/exx_lip.cpp
+++ b/source/module_ri/exx_lip.cpp
@@ -481,7 +481,7 @@ void Exx_Lip::b_cal( int ik, int iq, int ib)
 	}
 
 	std::complex<double> * const porter = new std::complex<double> [rho_basis->nrxx];
-	
+
 	for(size_t iw=0; iw< GlobalV::NLOCAL; ++iw)
 	{
 		const std::complex<double> * const phi_w = phi[iw];
@@ -495,7 +495,7 @@ void Exx_Lip::b_cal( int ik, int iq, int ib)
 		if( Conv_Coulomb_Pot_K::Ccp_Type::Ccp==info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf==info.ccp_type )
 			if((iq==iq_vecik) && (gzero_rank_in_pool==GlobalV::RANK_IN_POOL))							/// need to check while use k_point parallel
 				b0[iw] = b_w[rho_basis->ig_gge0];
-		
+
 		for( size_t ig=0; ig<rho_basis->npw; ++ig)
 			b_w[ig] *= recip_qkg2[ig];
 	}
@@ -634,12 +634,14 @@ void Exx_Lip::write_q_pack() const
 	if(!GlobalV::RANK_IN_POOL)
 	{
 		const std::string exx_q_pack = "exx_q_pack/";
-
+		int return_value=0;
 		const std::string command_mkdir = "test -d " + GlobalV::global_out_dir + exx_q_pack + " || mkdir " + GlobalV::global_out_dir + exx_q_pack;
-		system( command_mkdir.c_str() );	// Need to check
+        return_value = system(command_mkdir.c_str());
+        assert(return_value == 0);
 
-		const std::string command_kpoint = "test -f " + GlobalV::global_out_dir + exx_q_pack + GlobalV::global_kpoint_card + " || cp " + GlobalV::global_kpoint_card + " " + GlobalV::global_out_dir + exx_q_pack + GlobalV::global_kpoint_card;
-		system( command_kpoint.c_str() );	// Need to check
+        const std::string command_kpoint = "test -f " + GlobalV::global_out_dir + exx_q_pack + GlobalV::global_kpoint_card + " || cp " + GlobalV::global_kpoint_card + " " + GlobalV::global_out_dir + exx_q_pack + GlobalV::global_kpoint_card;
+        return_value = system(command_kpoint.c_str());
+		assert(return_value==0);
 
 		std::stringstream ss_wf_wg;
 		ss_wf_wg << GlobalV::global_out_dir << exx_q_pack << "wf_wg_" << GlobalV::MY_POOL;

From 05d462546207405973cef1583fdf4f8d960b2cf6 Mon Sep 17 00:00:00 2001
From: Denghui Lu <denghuilu@pku.edu.cn>
Date: Mon, 22 Jan 2024 10:21:02 +0800
Subject: [PATCH 07/17] output the default precision flag (#3496)

Co-authored-by: kirk0830 <67682086+kirk0830@users.noreply.github.com>
---
 source/module_io/write_input.cpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp
index b7e49901c8..256a7dce73 100644
--- a/source/module_io/write_input.cpp
+++ b/source/module_io/write_input.cpp
@@ -85,6 +85,7 @@ void Input::Print(const std::string &fn) const
     ModuleBase::GlobalFunc::OUTP(ofs, "cal_force", cal_force, "if calculate the force at the end of the electronic iteration");
     ModuleBase::GlobalFunc::OUTP(ofs, "out_freq_ion", out_freq_ion, "the frequency ( >= 0 ) of ionic step to output charge density and wavefunction. 0: output only when ion steps are finished");
     ModuleBase::GlobalFunc::OUTP(ofs, "device", device, "the computing device for ABACUS");
+    ModuleBase::GlobalFunc::OUTP(ofs, "precision", precision, "the computing precision for ABACUS");
 
     ofs << "\n#Parameters (2.PW)" << std::endl;
     ModuleBase::GlobalFunc::OUTP(ofs, "ecutwfc", ecutwfc, "#energy cutoff for wave functions");

From f91b2d61e0836c0794da6d58c0fdfd5d119d8071 Mon Sep 17 00:00:00 2001
From: YI Zeping <18586016708@163.com>
Date: Mon, 22 Jan 2024 03:56:44 +0000
Subject: [PATCH 08/17] Build: Improving CMake performance for finding LibXC
 and ELPA (#3478)

* Fix for finding LibXC and ELPA

* For compatibility to previous routines

* syntax fix for FindELPA.cmake

* Update cmake/FindELPA.cmake

Co-authored-by: Chun Cai <amoycaic@gmail.com>

* Using CMake interface as default for finding LibXC

* update docs

* fix for FindLibxc: changing imcompatible if statement

* fix for FindLibxc: changing imcompatible if statement

* fix for FindLibxc: changing imcompatible if statement

* update docs for installing pkg-config

* Update FindLibxc.cmake

* Update FindLibxc.cmake

* remove previous LibXC routine in CMakeLists.txt

Co-authored-by: Chun Cai <amoycaic@gmail.com>

* Update easy_install.md with Makefile-built LibXC supported

* Update easy_install.md to include different behavior in different version on finding ELPA

---------

Co-authored-by: Chun Cai <amoycaic@gmail.com>
---
 CMakeLists.txt                   |  7 ++---
 cmake/FindELPA.cmake             | 47 +++++++++++++++++++++-----------
 cmake/FindLibxc.cmake            | 36 ++++++++++++++++++++++++
 docs/quick_start/easy_install.md |  6 ++--
 4 files changed, 72 insertions(+), 24 deletions(-)
 create mode 100644 cmake/FindLibxc.cmake

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 99a6fc2650..0ee1e3bbc0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -536,11 +536,8 @@ if(DEFINED Libxc_DIR)
   set(ENABLE_LIBXC ON)
 endif()
 if(ENABLE_LIBXC)
-  find_package(Libxc REQUIRED HINTS
-    ${Libxc_DIR}/share/cmake/Libxc
-    ${Libxc_DIR}/lib/cmake/Libxc
-    ${Libxc_DIR}/lib64/cmake/Libxc
-  )
+  # use `cmake/FindLibxc.cmake` to detect Libxc installation with `pkg-config`
+  find_package(Libxc REQUIRED)
   message(STATUS "Found Libxc: version " ${Libxc_VERSION})
   if(${Libxc_VERSION} VERSION_LESS 5.1.7)
     message(FATAL_ERROR "LibXC >= 5.1.7 is required.")
diff --git a/cmake/FindELPA.cmake b/cmake/FindELPA.cmake
index 5769f7248c..328b75035f 100644
--- a/cmake/FindELPA.cmake
+++ b/cmake/FindELPA.cmake
@@ -7,34 +7,49 @@
 #  ELPA_INCLUDE_DIR - Where to find ELPA headers.
 #
 
-find_path(ELPA_INCLUDE_DIR
+find_package(PkgConfig)
+
+if(PKG_CONFIG_FOUND)
+  if(DEFINED ELPA_DIR)
+    string(APPEND CMAKE_PREFIX_PATH ";${ELPA_DIR}")
+  endif()
+  if(USE_OPENMP)
+    pkg_search_module(ELPA REQUIRED IMPORTED_TARGET GLOBAL elpa_openmp)
+  else()
+    pkg_search_module(ELPA REQUIRED IMPORTED_TARGET GLOBAL elpa)
+  endif()
+else()
+  find_path(ELPA_INCLUDE_DIRS
     elpa/elpa.h
     HINTS ${ELPA_DIR}
     PATH_SUFFIXES "include" "include/elpa"
     )
-if(USE_OPENMP)
-    find_library(ELPA_LIBRARY
-        NAMES elpa_openmp elpa
-        HINTS ${ELPA_DIR}
-        PATH_SUFFIXES "lib"
-        )
-else()
-    find_library(ELPA_LIBRARY
-        NAMES elpa
-        HINTS ${ELPA_DIR}
-        PATH_SUFFIXES "lib"
-        )
+  if(USE_OPENMP)
+    find_library(ELPA_LINK_LIBRARIES
+      NAMES elpa_openmp elpa
+      HINTS ${ELPA_DIR}
+      PATH_SUFFIXES "lib"
+      )
+  else()
+    find_library(ELPA_LINK_LIBRARIES
+      NAMES elpa
+      HINTS ${ELPA_DIR}
+      PATH_SUFFIXES "lib"
+      )
+  endif()
+  #message(
+  #  "ELPA : We need pkg-config to get all information about the elpa library")
 endif()
 
 # Handle the QUIET and REQUIRED arguments and
 # set ELPA_FOUND to TRUE if all variables are non-zero.
 include(FindPackageHandleStandardArgs)
-find_package_handle_standard_args(ELPA DEFAULT_MSG ELPA_LIBRARY ELPA_INCLUDE_DIR)
+find_package_handle_standard_args(ELPA DEFAULT_MSG ELPA_LINK_LIBRARIES ELPA_INCLUDE_DIRS)
 
 # Copy the results to the output variables and target.
 if(ELPA_FOUND)
-    set(ELPA_LIBRARIES ${ELPA_LIBRARY})
-    set(ELPA_INCLUDE_DIR ${ELPA_INCLUDE_DIR})
+    set(ELPA_LIBRARY ${ELPA_LINK_LIBRARIES})
+    set(ELPA_INCLUDE_DIR ${ELPA_INCLUDE_DIRS})
 
     if(NOT TARGET ELPA::ELPA)
         add_library(ELPA::ELPA UNKNOWN IMPORTED)
diff --git a/cmake/FindLibxc.cmake b/cmake/FindLibxc.cmake
new file mode 100644
index 0000000000..4a3c04cba7
--- /dev/null
+++ b/cmake/FindLibxc.cmake
@@ -0,0 +1,36 @@
+include(FindPackageHandleStandardArgs)
+
+if(DEFINED Libxc_DIR)
+  string(APPEND CMAKE_PREFIX_PATH ";${Libxc_DIR}")
+endif()
+# Using CMake interface as default.
+# NO REQUIRED here, otherwhile it would throw error
+# with no LibXC found.
+find_package(Libxc HINTS
+    ${Libxc_DIR}/share/cmake/Libxc
+    ${Libxc_DIR}/lib/cmake/Libxc
+    ${Libxc_DIR}/lib64/cmake/Libxc
+  )
+if(NOT TARGET Libxc::xc)
+  find_package(PkgConfig REQUIRED)
+  pkg_search_module(Libxc REQUIRED IMPORTED_TARGET GLOBAL libxc)
+  find_package_handle_standard_args(Libxc DEFAULT_MSG Libxc_LINK_LIBRARIES Libxc_INCLUDE_DIRS)
+endif()
+
+
+# Copy the results to the output variables and target.
+# if find_package() above works, Libxc::xc would be present and
+# below would be skipped.
+if(Libxc_FOUND AND NOT TARGET Libxc::xc)
+	set(Libxc_LIBRARY ${Libxc_LINK_LIBRARIES})
+	set(Libxc_LIBRARIES ${Libxc_LIBRARY})
+	set(Libxc_INCLUDE_DIR ${Libxc_INCLUDE_DIRS})
+	add_library(Libxc::xc UNKNOWN IMPORTED)
+	set_target_properties(Libxc::xc PROPERTIES
+		IMPORTED_LOCATION "${Libxc_LIBRARY}"
+		INTERFACE_INCLUDE_DIRECTORIES "${Libxc_INCLUDE_DIR}")
+endif()
+
+set(CMAKE_REQUIRED_INCLUDES ${CMAKE_REQUIRED_INCLUDES} ${Libxc_INCLUDE_DIR})
+
+mark_as_advanced(Libxc_INCLUDE_DIR Libxc_LIBRARY)
diff --git a/docs/quick_start/easy_install.md b/docs/quick_start/easy_install.md
index 32ae5b87d3..f2173d569e 100644
--- a/docs/quick_start/easy_install.md
+++ b/docs/quick_start/easy_install.md
@@ -28,7 +28,7 @@ These requirements support the calculation of plane-wave basis in ABACUS. For LC
 Some of these packages can be installed with popular package management system, such as `apt` and `yum`:
 
 ```bash
-sudo apt update && sudo apt install -y libopenblas-openmp-dev liblapack-dev libscalapack-mpi-dev libelpa-dev libfftw3-dev libcereal-dev libxc-dev g++ make cmake bc git
+sudo apt update && sudo apt install -y libopenblas-openmp-dev liblapack-dev libscalapack-mpi-dev libelpa-dev libfftw3-dev libcereal-dev libxc-dev g++ make cmake bc git pkgconf
 ```
 
 > Installing ELPA by apt only matches requirements on Ubuntu 22.04. For earlier linux distributions, you should build ELPA from source.
@@ -111,12 +111,12 @@ Here, 'build' is the path for building ABACUS; and '-D' is used for setting up s
   - `LAPACK_DIR`: Path to OpenBLAS library `libopenblas.so`(including BLAS and LAPACK)
   - `SCALAPACK_DIR`: Path to ScaLAPACK library `libscalapack.so`
   - `ELPA_DIR`: Path to ELPA install directory; should be the folder containing 'include' and 'lib'.
-  > Note: If you install ELPA from source, please add a symlink to avoid the additional include file folder with version name: `ln -s elpa/include/elpa-2021.05.002/elpa elpa/include/elpa`. This is a known behavior of ELPA.
+  > Note: In ABACUS v3.5.1 or earlier, if you install ELPA from source , please add a symlink to avoid the additional include file folder with version name: `ln -s elpa/include/elpa-2021.05.002/elpa elpa/include/elpa` to help the build system find ELPA headers.
 
   - `FFTW3_DIR`: Path to FFTW3.
   - `CEREAL_INCLUDE_DIR`: Path to the parent folder of `cereal/cereal.hpp`. Will download from GitHub if absent.
   - `Libxc_DIR`: (Optional) Path to Libxc.
-  > Note: Building Libxc from source with Makefile does NOT support using it in CMake here. Please compile Libxc with CMake instead.
+  > Note: In ABACUS v3.5.1 or earlier, Libxc built from source with Makefile is NOT supported; please compile Libxc with CMake instead.
   - `LIBRI_DIR`: (Optional) Path to LibRI.
   - `LIBCOMM_DIR`: (Optional) Path to LibComm.
 

From 2189ba160ab8c0620d78b044d304a1941381be07 Mon Sep 17 00:00:00 2001
From: wqzhou <33364058+WHUweiqingzhou@users.noreply.github.com>
Date: Mon, 22 Jan 2024 14:58:08 +0800
Subject: [PATCH 09/17] Docs: correct some docs about mp2 smearing method
 (#3533)

* correct some docs about mp2 smearing method

* add docs about mv method
---
 docs/advanced/input_files/input-main.md |   2 +
 source/module_elecstate/occupy.cpp      | 414 +-----------------------
 2 files changed, 8 insertions(+), 408 deletions(-)

diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
index fb721a7486..1082f4a88b 100644
--- a/docs/advanced/input_files/input-main.md
+++ b/docs/advanced/input_files/input-main.md
@@ -945,6 +945,8 @@ calculations.
   - **fixed**: fixed occupations (available for non-coductors only)
   - **gauss** or **gaussian**: Gaussian smearing method.
   - **mp**: methfessel-paxton smearing method; recommended for metals.
+  - **mp2**: 2-nd methfessel-paxton smearing method; recommended for metals.
+  - **mv** or **cold**: marzari-vanderbilt smearing method.
   - **fd**: Fermi-Dirac smearing method: $f=1/\{1+\exp[(E-\mu)/kT]\}$ and smearing_sigma below is the temperature $T$ (in Ry).
 - **Default**: gauss
 
diff --git a/source/module_elecstate/occupy.cpp b/source/module_elecstate/occupy.cpp
index 80918dd3f1..e896aae4e8 100644
--- a/source/module_elecstate/occupy.cpp
+++ b/source/module_elecstate/occupy.cpp
@@ -79,6 +79,12 @@ void Occupy::decision(const std::string &name, const std::string &smearing_metho
         {
             gaussian_type = 2; // 2nd Methfessel-Paxton method.
         }
+        else if (smearing_method == "mp3")
+        {
+            // acually any order Methfessel-Paxton method can be supported in Occupy::w1gauss()
+            // however the parameter is string instead of int
+            ModuleBase::WARNING_QUIT("occupy", "Some refactor of smearing shoule be done before supporting any order of Methfessel-Paxton method!");
+        }
 
         else if (smearing_method == "marzari-vanderbilt" || smearing_method == "cold" || smearing_method == "mv")
         {
@@ -597,411 +603,3 @@ double Occupy::w1gauss(const double &x, const int n)
 
     return w1;
 } // end function w1gauss
-
-/*
-void Occupy::tweights(const int nks,const int nspin,const int nband,const double &nelec,
-                      const int ntetra,const ModuleBase::matrix &tetra, double **ekb, double &ef, ModuleBase::matrix
-&wg)
-{
-    //===================================================================
-    // calculates weights with the tetrahedron method (Bloechl version)
-    // integer :: nks, nspin, GlobalV::NBANDS, ntetra, tetra (4, ntetra)
-    //===================================================================
-
-    double e1, e2, e3, e4, c1, c2, c3, c4, dosef;
-    int ik, ibnd, nt, nk, ns, i, kp1, kp2, kp3, kp4;
-
-    double etetra[4];
-    int itetra[4];
-
-    // Calculate the Fermi energy ef
-    efermit(ekb, GlobalV::NBANDS, nks, nelec, nspin, ntetra, tetra, ef);
-
-    for (ik = 0;ik < nks;ik++)
-    {
-        for (ibnd = 0;ibnd < nband;ibnd++)
-        {
-            wg(ik, ibnd) = 0.0;
-        } // enddo
-    } // enddo
-
-    for (ns = 0;ns < nspin;ns++)
-    {
-        //==================================================================
-        // nk is used to select k-points with up (ns=1) or down (ns=2) spin
-        //==================================================================
-        if (ns == 1)
-        {
-            nk = 0;
-        }
-        else
-        {
-            nk = nks / 2;
-        }
-
-        for (nt = 0;nt < ntetra;nt++)
-        {
-            for (ibnd = 0;ibnd < GlobalV::NBANDS;ibnd++)
-            {
-                //======================================================
-                // etetra are the energies at the vertexes of the nt-th
-                // tetrahedron
-                //======================================================
-                for (i = 0;i < 4;i++)
-                {
-                    etetra [i] = ekb[static_cast<int>( tetra(nt,i) ) + nk][ibnd];
-                }
-
-                itetra[0] = 0;
-
-                ModuleBase::hpsort(4, etetra, itetra);
-
-                //===============================================
-                // ...sort in ascending order: e1 < e2 < e3 < e4
-                //===============================================
-                e1 = etetra [0];
-                e2 = etetra [1];
-                e3 = etetra [2];
-                e4 = etetra [3];
-
-                //==============================================================
-                // kp1-kp4 are the irreducible k-points corresponding to e1-e4
-                //==============================================================
-
-                kp1 = static_cast<int>( tetra(nt,itetra[0]) )+ nk;
-                kp2 = static_cast<int>( tetra(nt,itetra[1]) )+ nk;
-                kp3 = static_cast<int>( tetra(nt,itetra[2]) )+ nk;
-                kp4 = static_cast<int>( tetra(nt,itetra[3]) )+ nk;
-
-                //======================
-                // calculate weights wg
-                //======================
-                if (ef >= e4)
-                {
-                    wg(kp1, ibnd) = wg(kp1, ibnd) + 0.250 / ntetra;
-                    wg(kp2, ibnd) = wg(kp2, ibnd) + 0.250 / ntetra;
-                    wg(kp3, ibnd) = wg(kp3, ibnd) + 0.250 / ntetra;
-                    wg(kp4, ibnd) = wg(kp4, ibnd) + 0.250 / ntetra;
-                }
-                else if (ef < e4 && ef >= e3)
-                {
-                    c4 = 0.250 / ntetra * pow(e4 - ef, 3) / (e4 - e1) / (e4 - e2)
-                         / (e4 - e3);
-                    dosef = 3.0 / ntetra * (e4 - ef) * (e4 - ef) / (e4 - e1) / (e4 - e2)
-                            / (e4 - e3);
-                    wg(kp1, ibnd) = wg(kp1, ibnd) + 0.250 / ntetra - c4 *
-                                    (e4 - ef) / (e4 - e1) + dosef * (e1 + e2 + e3 + e4 - 4.0 * ekb[kp1][ibnd]) / 40.0;
-                    wg(kp2, ibnd) = wg(kp2, ibnd) + 0.250 / ntetra - c4 *
-                                    (e4 - ef) / (e4 - e2) + dosef * (e1 + e2 + e3 + e4 - 4.0 * ekb[kp2][ibnd]) / 40.0;
-                    wg(kp3, ibnd) = wg(kp3, ibnd) + 0.250 / ntetra - c4 *
-                                    (e4 - ef) / (e4 - e3) + dosef * (e1 + e2 + e3 + e4 - 4.0 * ekb[kp3][ibnd]) / 40.0;
-                    wg(kp4, ibnd) = wg(kp4, ibnd) + 0.250 / ntetra - c4 *
-                                    (4.0 - (e4 - ef) * (1.0 / (e4 - e1) + 1.0 / (e4 - e2)
-                                                        + 1.0 / (e4 - e3))) + dosef * (e1 + e2 + e3 + e4 - 4.0 *
-                                                                                       ekb[kp4][ibnd]) / 40.0;
-                }
-
-                else if (ef < e3 && ef >= e2)
-                {
-                    c1 = 0.250 / ntetra * (ef - e1) * (ef - e1) / (e4 - e1) / (e3 - e1);
-                    c2 = 0.250 / ntetra * (ef - e1) * (ef - e2) * (e3 - ef)
-                         / (e4 - e1) / (e3 - e2) / (e3 - e1);
-                    c3 = 0.250 / ntetra * (ef - e2) * (ef - e2) * (e4 - ef) / (e4 - e2)
-                         / (e3 - e2) / (e4 - e1);
-                    dosef = 1.0 / ntetra / (e3 - e1) / (e4 - e1) * (3.0 *
-                            (e2 - e1) + 6.0 * (ef - e2) - 3.0 * (e3 - e1 + e4 - e2)
-                            * (ef - e2) * (ef - e2) / (e3 - e2) / (e4 - e2));
-                    wg(kp1, ibnd) = wg(kp1, ibnd) + c1 + (c1 + c2) * (e3 - ef)
-                                    / (e3 - e1) + (c1 + c2 + c3) * (e4 - ef) / (e4 - e1) + dosef *
-                                    (e1 + e2 + e3 + e4 - 4.0 * ekb[kp1][ibnd]) / 40.0;
-                    wg(kp2, ibnd) = wg(kp2, ibnd) + c1 + c2 + c3 + (c2 + c3)
-                                    * (e3 - ef) / (e3 - e2) + c3 * (e4 - ef) / (e4 - e2) + dosef *
-                                    (e1 + e2 + e3 + e4 - 4.0 * ekb[kp2][ibnd]) / 40.0;
-                    wg(kp3, ibnd) = wg(kp3, ibnd) + (c1 + c2) * (ef - e1)
-                                    / (e3 - e1) + (c2 + c3) * (ef - e2) / (e3 - e2) + dosef *
-                                    (e1 + e2 + e3 + e4 - 4.0 * ekb[kp3][ibnd]) / 40.0;
-                    wg(kp4, ibnd) = wg(kp4, ibnd) + (c1 + c2 + c3) * (ef - e1)
-                                    / (e4 - e1) + c3 * (ef - e2) / (e4 - e2) + dosef * (e1 + e2 +
-                                            e3 + e4 - 4.0 * ekb[kp4][ibnd]) / 40.0;
-                }
-                else if (ef < e2 && ef >= e1)
-                {
-                    c4 = 0.250 / ntetra * (ef - e1) * (ef - e1) * (ef - e1) / (e2 - e1) /
-                         (e3 - e1) / (e4 - e1);
-                    dosef = 3.0 / ntetra * (ef - e1) * (ef - e1) / (e2 - e1) / (e3 - e1)
-                            / (e4 - e1);
-                    wg(kp1, ibnd) = wg(kp1, ibnd) + c4 * (4.0 - (ef - e1)
-                                                          * (1.0 / (e2 - e1) + 1.0 / (e3 - e1) + 1.0 / (e4 - e1)))
-                                    + dosef * (e1 + e2 + e3 + e4 - 4.0 * ekb[kp1][ibnd]) / 40.0;
-                    wg(kp2, ibnd) = wg(kp2, ibnd) + c4 * (ef - e1) / (e2 - e1)
-                                    + dosef * (e1 + e2 + e3 + e4 - 4.0 * ekb[kp2][ibnd]) / 40.0;
-                    wg(kp3, ibnd) = wg(kp3, ibnd) + c4 * (ef - e1) / (e3 - e1)
-                                    + dosef * (e1 + e2 + e3 + e4 - 4.0 * ekb[kp3][ibnd]) / 40.0;
-                    wg(kp4, ibnd) = wg(kp4, ibnd) + c4 * (ef - e1) / (e4 - e1)
-                                    + dosef * (e1 + e2 + e3 + e4 - 4.0 * ekb[kp4][ibnd]) / 40.0;
-                } // endif
-            } // enddo
-        } // enddo
-    } // enddo
-
-    //=====================================================================
-    // add correct spin normalization : 2 for LDA, 1 for LSDA calculations
-    //=====================================================================
-    for (ik = 0;ik < nks;ik++)
-    {
-        for (ibnd = 0;ibnd < GlobalV::NBANDS;ibnd++)
-        {
-            wg(ik, ibnd) = wg(ik, ibnd) * 2.0 / nspin;
-        }
-    }
-    return;
-} // end subroutine tweights
-*/
-
-/*
-double Occupy::wsweight(const ModuleBase::Vector3<double> &r, ModuleBase::Vector3<double> *rws,const int nrws)
-{
-    //============================================================
-    // integer ir, nreq, nrws
-    // real(kind=dp) r(3), rrt, ck, eps, rws(0:3,nrws), wsweight
-    // parameter (eps=1.0e-6)
-    //============================================================
-    const double eps = 1.0e-6;
-
-    int nreq = 1;
-
-    for (int ir = 0;ir < nrws;ir++)
-    {
-        const double rrt = r * rws[ir];
-        const double ck = rrt - rws[ir].x;
-        //	rrt = r[1]*rws(1,ir) + r[2]*rws(2,ir) + r[3]*rws(3,ir);
-        //	ck = rrt-rws(0,ir);
-
-        if (ck > eps)
-        {
-            break;
-        }
-
-        if (std::abs(ck) < eps)
-        {
-            nreq++;
-        }
-    } // end do
-
-    const double wswe = 1.0 / nreq;
-
-    return wswe;
-} // end function wsweight
-*/
-
-/*
-void Occupy::efermit(double** ekb,const int nband,const int nks,const double &nelec,const int nspin,
-                     const int ntetra,const ModuleBase::matrix &tetra, double &ef)
-{
-    //=======================================================
-    // Finds the Fermi energy - tetrahedron method (Bloechl)
-    // the transformation Ry to eV
-    //=======================================================
-
-    // parameter :
-    const int maxiter = 300;
-    const double eps = 1.0e-10;
-
-    double efbetter;
-
-    //===================================
-    // nlw : the minimum energy band
-    // elw : the lower limit of the fermi ener
-    // eup : the upper limit of the fermi ener
-    // external sumkt
-    // find bounds for the Fermi energy.
-    //===================================
-    const int nlw = max(  1, static_cast<int>( (nelec / 2.0 - 5.0) )  );
-    double elw = ekb[nlw][0];
-    double eup = ekb[0][GlobalV::NBANDS-1];
-
-    for (int ik = 1;ik < nks;ik++)// do ik = 2, nks
-    {
-        elw = min(elw, ekb[ik][nlw]);
-        eup = max(eup, ekb[ik][GlobalV::NBANDS-1]);
-    }
-    for (int ik = 1;ik < nks;ik++)// do ik = 2, nks
-    {
-        elw = min(elw, ekb[ik][nlw]);
-        eup = max(eup, ekb[ik][GlobalV::NBANDS-1]);
-    }
-
-    //===============================
-    // Bisection method
-    // the number of states with eup
-    // the number of states with elw
-    //===============================
-    const double sumkup = sumkt(ekb, GlobalV::NBANDS, nks, nspin, ntetra, tetra, eup);
-    const double sumklw = sumkt(ekb, GlobalV::NBANDS, nks, nspin, ntetra, tetra, elw);
-
-    GlobalV::ofs_running << "\n sumkup = " << sumkup;
-    GlobalV::ofs_running << "\n sumklw = " << sumklw << std::endl;
-
-    if ((sumkup - nelec) < - eps || (sumklw - nelec) > eps)
-    {
-        ModuleBase::WARNING("efermit","unexpected error.");
-    }
-
-    double better = 1.0e+10;
-
-    bool converge = false;
-
-    double sumkmid = 0.0;
-    for (int iter = 0;iter < maxiter;iter++)
-    {
-        // the number of states with ef
-        ef = (eup + elw) / 2.0;
-        sumkmid = sumkt(ekb, GlobalV::NBANDS, nks, nspin, ntetra, tetra, ef);
-
-        if (std::abs(sumkmid - nelec) < better)
-        {
-            better = std::abs(sumkmid - nelec);
-            efbetter = ef;
-        }
-
-        // converged
-        if (std::abs(sumkmid - nelec) < eps)
-        {
-            converge = true;
-            break;
-        }
-        else if ((sumkmid - nelec) < - eps)
-        {
-            elw = ef;
-        }
-        else
-        {
-            eup = ef;
-        }
-    }
-    if (!converge)
-    {
-        // unconverged exit:
-        // the best available ef is used . Needed in some difficult cases
-        ef = efbetter;
-        sumkmid = sumkt(ekb, GlobalV::NBANDS, nks, nspin, ntetra, tetra, ef);
-    }
-
-    //==============================================================
-    // Check if Fermi level is above any of the highest eigenvalues
-    //==============================================================
-    for (int ik = 0;ik < nks;ik++)
-    {
-        if (ef > ekb[ik][GlobalV::NBANDS-1] + 1.e-4)
-        {
-            std::cout << "\n ef = " << ef;
-        }
-    }
-    return;
-} // end subroutine efermit
-*/
-
-/*
-double Occupy::sumkt(double** ekb,const int nband,const int nks,const int nspin,const int ntetra,
-                     const ModuleBase::matrix &tetra,const double &e)
-{
-    double etetra[4];
-    double sum = 0.0;
-
-    int nk = 0 ;
-    for (int ns = 0; ns < nspin;ns++)
-    {
-        //==================================================================
-        // nk is used to select k-points with up (ns=1) or down (ns=2) spin
-        //==================================================================
-        if (ns == 1)
-        {
-            nk = 0;
-        }
-        else
-        {
-            nk = nks / 2;
-        }
-
-        for (int nt = 0; nt < ntetra; nt++)
-        {
-            for (int ibnd = 0; ibnd < GlobalV::NBANDS; ibnd++)
-            {
-                //======================================================
-                // etetra are the energies at the vertexes of the nt-th
-                // tetrahedron
-                //======================================================
-                for (int i = 0; i < 4; i++)
-                {
-                    etetra [i] = ekb[  static_cast<int>( (tetra(i, nt) + nk) )][ ibnd  ];
-                }
-
-                piksort(4, etetra);
-                //===========================================
-                //sort in ascending order: e1 < e2 < e3 < e4
-                //===========================================
-                const double e1 = etetra [0];
-                const double e2 = etetra [1];
-                const double e3 = etetra [2];
-                const double e4 = etetra [3];
-
-                //===============================================
-                // calculate sum over k of the integrated charge
-                //===============================================
-                if (e >= e4)
-                {
-                    sum += 1.0 / ntetra;
-                }
-                else if (e < e4 && e >= e3)
-                {
-                    sum += 1.0 / ntetra * (1.0 - pow((e4 - e), 3) / (e4 - e1)
-                                           / (e4 - e2) / (e4 - e3));
-                }
-                else if (e < e3 && e >= e2)
-                {
-                    sum += 1.0 / ntetra / (e3 - e1) / (e4 - e1) *
-                           ((e2 - e1) * (e2 - e1) + 3.0 * (e2 - e1) * (e - e2) +
-                            3.0 * (e - e2) * (e - e2) - (e3 - e1 + e4 - e2) /
-                            (e3 - e2) / (e4 - e2) * pow((e - e2), 3));
-                }
-                else if (e < e2 && e >= e1)
-                {
-                    sum += 1.0 / ntetra * pow((e - e1), 3) /
-                           (e2 - e1) / (e3 - e1) / (e4 - e1);
-                }
-            }//ibnd
-        }//nt
-    }//ns
-
-// add correct spin normalization : 2 for LDA, 1 for LSDA calculations
-    sum *= 2.0 / nspin;
-    return sum;
-} // end function sumkt
-*/
-
-/*
-void Occupy::piksort(const int n, double *a)
-{
-    int i;
-    bool b = true;
-    for (int j = 1;j < n;j++) // do j = 2, n
-    {
-        const double temp = a [j];
-        for (i = j - 1;i >= 0;i--)  // do i = j - 1, 1, - 1
-        {
-            if (a [i] <= temp)
-            {
-                b = false;
-                break;
-            }
-            a [i + 1] = a [i];
-        }
-        if (b)
-        {
-            i = 0;
-        }
-        a [i + 1] = temp;
-    }
-    return;
-} //end subroutine piksort
-*/

From 94525b56a7083a5b03960b674dfc2a55383bf545 Mon Sep 17 00:00:00 2001
From: Wenfei Li <38569667+wenfei-li@users.noreply.github.com>
Date: Mon, 22 Jan 2024 15:00:45 +0800
Subject: [PATCH 10/17] Feature : printing band density (#3501)

Co-authored-by: wenfei-li <liwenfei@gmail.com>
Co-authored-by: wqzhou <33364058+WHUweiqingzhou@users.noreply.github.com>
---
 source/module_esolver/esolver_ks_pw.cpp | 52 +++++++++++++++++++++++++
 source/module_io/input.cpp              | 44 +++++++++++++++++++++
 source/module_io/input.h                |  2 +
 3 files changed, 98 insertions(+)

diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
index 2e81938ae2..0bdf75045f 100644
--- a/source/module_esolver/esolver_ks_pw.cpp
+++ b/source/module_esolver/esolver_ks_pw.cpp
@@ -894,6 +894,58 @@ void ESolver_KS_PW<T, Device>::afterscf(const int istep)
                             this->kspw_psi[0].get_pointer() - this->kspw_psi[0].get_psi_bias(),
                             this->psi[0].size());
     }
+
+    if(INPUT.band_print_num > 0)
+    {
+        std::complex<double> * wfcr = new std::complex<double>[this->pw_rho->nxyz];
+        double * rho_band = new double [this->pw_rho->nxyz];
+        for(int i = 0; i < this->pw_rho->nxyz; i++)
+        {
+            rho_band[i] = 0.0;
+        }
+
+        for(int i = 0; i < INPUT.band_print_num; i++)
+        {
+            int ib = INPUT.bands_to_print[i];
+            for(int ik = 0; ik < this->kv.nks; ik++)
+            {
+                this->psi->fix_k(ik);
+                this->pw_wfc->recip_to_real(this->ctx,&psi[0](ib,0),wfcr,ik);
+
+                double w1 = static_cast<double>(this->kv.wk[ik] / GlobalC::ucell.omega);
+
+                for(int i = 0; i < this->pw_rho->nxyz; i++)
+                {
+                    rho_band[i] += std::norm(wfcr[i]) * w1;
+                }
+            }
+
+            std::stringstream ssc;
+            ssc << GlobalV::global_out_dir << "band" << ib << ".cube";     
+
+            ModuleIO::write_rho
+            (
+#ifdef __MPI
+                this->pw_big->bz,
+                this->pw_big->nbz,
+                this->pw_big->nplane,
+                this->pw_big->startz_current,
+#endif
+                rho_band,
+                0,
+                GlobalV::NSPIN,
+                0,
+                ssc.str(),
+                this->pw_rho->nx,
+                this->pw_rho->ny,
+                this->pw_rho->nz,
+                0.0,
+                &(GlobalC::ucell),
+                11);
+        }
+        delete[] wfcr;
+        delete[] rho_band;
+    }
 }
 
 template <typename T, typename Device>
diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp
index abe8f28ebd..7a8d00a838 100644
--- a/source/module_io/input.cpp
+++ b/source/module_io/input.cpp
@@ -327,6 +327,8 @@ void Input::Default(void)
 
     out_bandgap = 0; // QO added for bandgap printing
 
+    band_print_num = 0;
+
     deepks_out_labels = 0; // caoyu added 2020-11-24, mohan added 2021-01-03
     deepks_scf = 0;
     deepks_bandgap = 0;
@@ -1327,6 +1329,14 @@ bool Input::Read(const std::string& fn)
         {
             read_bool(ifs, out_chg);
         }
+        else if (strcmp("band_print_num", word) == 0)
+        {
+            read_value(ifs, band_print_num);
+        }
+        else if (strcmp("bands_to_print", word) == 0)
+        {
+            ifs.ignore(150, '\n');
+        }
         else if (strcmp("out_dm", word) == 0)
         {
             read_bool(ifs, out_dm);
@@ -2369,6 +2379,29 @@ bool Input::Read(const std::string& fn)
         ModuleBase::WARNING_QUIT("Input", "The ntype in INPUT is not equal to the ntype counted in STRU, check it.");
     }
 
+    if(band_print_num > 0)
+    {
+        bands_to_print.resize(band_print_num);
+        ifs.clear();
+        ifs.seekg(0); // move to the beginning of the file
+        ifs.rdstate();
+        while (ifs.good())
+        {
+            ifs >> word1;
+            if (ifs.eof() != 0)
+                break;
+            strtolower(word1, word); // convert uppercase std::string to lower case; word1 --> word
+
+            if (strcmp("bands_to_print", word) == 0)
+            {
+                for(int i = 0; i < band_print_num; i ++)
+                {
+                    ifs >> bands_to_print[i];
+                }
+            }
+        }
+    }
+
     //----------------------------------------------------------
     //       DFT+U    Xin Qu  added on 2020-10-29
     //----------------------------------------------------------
@@ -3524,6 +3557,17 @@ void Input::Bcast()
     Parallel_Common::bcast_bool(restart_save);  // Peize Lin add 2020.04.04
     Parallel_Common::bcast_bool(restart_load);  // Peize Lin add 2020.04.04
 
+    Parallel_Common::bcast_int(band_print_num);
+    if(GlobalV::MY_RANK != 0)
+    {
+        bands_to_print.resize(band_print_num);
+    }
+
+    for(int i = 0; i < band_print_num; i++)
+    {
+        Parallel_Common::bcast_int(bands_to_print[i]);
+    }
+
     //-----------------------------------------------------------------------------------
     // DFT+U (added by Quxin 2020-10-29)
     //-----------------------------------------------------------------------------------
diff --git a/source/module_io/input.h b/source/module_io/input.h
index 6393483cb4..d0344b8944 100644
--- a/source/module_io/input.h
+++ b/source/module_io/input.h
@@ -259,6 +259,8 @@ class Input
     bool out_chg; // output charge density. 0: no; 1: yes
     bool out_dm; // output density matrix.
     bool out_dm1;
+    int band_print_num;
+    std::vector<int> bands_to_print;
     int out_pot; // yes or no
     int out_wfc_pw; // 0: no; 1: txt; 2: dat
     bool out_wfc_r; // 0: no; 1: yes

From ebbb6f3e47c28d8568fa5e8184f8b978b9e201be Mon Sep 17 00:00:00 2001
From: wqzhou <33364058+WHUweiqingzhou@users.noreply.github.com>
Date: Tue, 23 Jan 2024 15:31:00 +0800
Subject: [PATCH 11/17] add some docs for PR#3501 (#3537)

---
 docs/advanced/input_files/input-main.md | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
index 1082f4a88b..13a555b23e 100644
--- a/docs/advanced/input_files/input-main.md
+++ b/docs/advanced/input_files/input-main.md
@@ -146,6 +146,8 @@
     - [out\_ndigits](#out_ndigits)
     - [out\_interval](#out_interval)
     - [out\_element\_info](#out_element_info)
+    - [band\_print\_num](#band_print_num)
+    - [bands\_to\_print](#bands_to_print)
     - [restart\_save](#restart_save)
     - [restart\_load](#restart_load)
     - [rpa](#rpa)
@@ -1601,6 +1603,20 @@ These variables are used to control the output of properties.
 - **Description**: Control the interval for printing Mulliken population analysis, $r(R)$, $H(R)$, $S(R)$, $T(R)$, $dH(R)$, $H(k)$, $S(k)$ and $wfc(k)$ matrices during molecular dynamics calculations. Check input parameters [out_mul](#out_mul), [out_mat_r](#out_mat_r), [out_mat_hs2](#out_mat_hs2), [out_mat_t](#out_mat_t), [out_mat_dh](#out_mat_dh), [out_mat_hs](#out_mat_hs) and [out_wfc_lcao](#out_wfc_lcao) for more information, respectively.
 - **Default**: 1
 
+### band_print_num
+
+- **Type**: Integer
+- **Availability**: PW basis
+- **Description**: If you want to plot a partial charge density contributed from some chosen bands. `band_print_num` define the number of band list. The result can be found in "band*.cube".
+- **Default**: 0
+
+### bands_to_print
+
+- **Type**: vector
+- **Availability**: band_print_num > 0
+- **Description**: define which band you want to choose for partial charge density.
+- **Default**: []
+
 ### out_element_info
 
 - **Type**: Boolean

From 30f2de895620aa2b527e175f8ba84659bcef8a11 Mon Sep 17 00:00:00 2001
From: wqzhou <33364058+WHUweiqingzhou@users.noreply.github.com>
Date: Tue, 23 Jan 2024 17:38:34 +0800
Subject: [PATCH 12/17] Feature: enable restart charge density mixing during
 SCF (#3542)

* add a new parameter mixing_restart

* do not update rho if iter==mixing_restart

* do not update rho if iter==mixing_restart-1

* reset mix and rho_mdata if iter==mixing_restart

* fix SCF exit directly since drho=0 if iter=GlobalV::MIXING_RESTART

* re-set_mixing in eachiterinit for PW and LCAO

* enable SCF restarts in esolver_ks::RUN

* add some UnitTests

* add some Docs
---
 docs/advanced/input_files/input-main.md    |  7 +++++++
 source/module_base/global_variable.cpp     |  1 +
 source/module_base/global_variable.h       |  1 +
 source/module_esolver/esolver_ks.cpp       | 19 ++++++++++++++++---
 source/module_esolver/esolver_ks_lcao.cpp  | 13 ++++++++++++-
 source/module_esolver/esolver_ks_pw.cpp    | 14 ++++++++++++--
 source/module_io/input.cpp                 |  6 ++++++
 source/module_io/input.h                   |  1 +
 source/module_io/input_conv.cpp            |  1 +
 source/module_io/parameter_pool.cpp        |  4 ++++
 source/module_io/test/input_conv_test.cpp  |  1 +
 source/module_io/test/input_test_para.cpp  |  1 +
 source/module_io/test/write_input_test.cpp | 11 +++++++----
 source/module_io/write_input.cpp           |  1 +
 14 files changed, 71 insertions(+), 10 deletions(-)

diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
index 13a555b23e..42a4b3638a 100644
--- a/docs/advanced/input_files/input-main.md
+++ b/docs/advanced/input_files/input-main.md
@@ -1005,6 +1005,13 @@ We recommend the following options:
   For systems that are difficult to converge, one could try increasing the value of 'mixing_ndim' to enhance the stability of the self-consistent field (SCF) calculation.
 - **Default**: 8
 
+### mixing_restart
+
+- **Type**: Integer
+- **Description**: At `mixing_restart`-th iteration, SCF will restart by using output charge density from perivos iteration as input charge density directly, and start a new mixing. `mixing_restart=0|1` means SCF starts from scratch.
+  
+- **Default**: 0
+
 ### mixing_gg0
 
 - **Type**: Real
diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp
index 696bcd6088..eb0dc636e2 100644
--- a/source/module_base/global_variable.cpp
+++ b/source/module_base/global_variable.cpp
@@ -248,6 +248,7 @@ std::string of_kernel_file = "WTkernel.txt";
 std::string MIXING_MODE = "broyden";
 double MIXING_BETA = 0.7;
 int MIXING_NDIM = 8;
+int MIXING_RESTART = 0;
 double MIXING_GG0 = 1.00;
 double MIXING_BETA_MAG = 1.6;
 double MIXING_GG0_MAG = 1.00;
diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h
index b1fbb1748d..9fa7dc1c8a 100644
--- a/source/module_base/global_variable.h
+++ b/source/module_base/global_variable.h
@@ -277,6 +277,7 @@ extern std::string of_kernel_file; // The name of WT kernel file.
 extern std::string MIXING_MODE;
 extern double MIXING_BETA;
 extern int MIXING_NDIM;
+extern int MIXING_RESTART;
 extern double MIXING_GG0;
 extern bool MIXING_TAU;
 extern double MIXING_BETA_MAG;
diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp
index d9877a0969..b4fcd83375 100644
--- a/source/module_esolver/esolver_ks.cpp
+++ b/source/module_esolver/esolver_ks.cpp
@@ -410,7 +410,7 @@ namespace ModuleESolver
                         }
                     }
 
-                    this->conv_elec = (drho < this->scf_thr);
+                    this->conv_elec = (drho < this->scf_thr && iter!=GlobalV::MIXING_RESTART);
 
                     // If drho < hsolver_error in the first iter or drho < scf_thr, we do not change rho.
                     if (drho < hsolver_error || this->conv_elec)
@@ -436,8 +436,16 @@ namespace ModuleESolver
                         //     }
                         //     p_chgmix->auto_set(bandgap_for_autoset, GlobalC::ucell);
                         // }
-                        
-                        p_chgmix->mix_rho(pelec->charge);
+                        // mixing will restart after GlobalV::MIXING_RESTART steps
+                        // So, GlobalV::MIXING_RESTART=1 means mix from scratch
+                        if (GlobalV::MIXING_RESTART > 0 && iter == GlobalV::MIXING_RESTART - 1)
+                        {
+                            // do not mix charge density
+                        }
+                        else
+                        {
+                            p_chgmix->mix_rho(pelec->charge); // update chr->rho by mixing
+                        }
                         if (GlobalV::SCF_THR_TYPE == 2) pelec->charge->renormalize_rho(); // renormalize rho in R-space would induce a error in K-space
                         //----------charge mixing done-----------
                     }
@@ -467,6 +475,11 @@ namespace ModuleESolver
                     bool stop = this->do_after_converge(iter);
                     if(stop) break;
                 }
+                // notice for restart
+                if (GlobalV::MIXING_RESTART > 0 && iter == GlobalV::MIXING_RESTART - 1)
+                {
+                    std::cout<<"SCF restart after this step!"<<std::endl;
+                }
             }
             afterscf(istep);
 
diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
index a7d0a5abdd..4c4d6c342f 100644
--- a/source/module_esolver/esolver_ks_lcao.cpp
+++ b/source/module_esolver/esolver_ks_lcao.cpp
@@ -491,8 +491,19 @@ namespace ModuleESolver
     template <typename TK, typename TR>
     void ESolver_KS_LCAO<TK, TR>::eachiterinit(const int istep, const int iter)
 {
-    if (iter == 1)
+    if (iter == 1 || iter == GlobalV::MIXING_RESTART)
+    {
+        if (iter == GlobalV::MIXING_RESTART) // delete mixing and re-construct it to restart 
+        {
+            this->p_chgmix->set_mixing(GlobalV::MIXING_MODE,
+                                GlobalV::MIXING_BETA,
+                                GlobalV::MIXING_NDIM,
+                                GlobalV::MIXING_GG0,
+                                GlobalV::MIXING_TAU,
+                                GlobalV::MIXING_BETA_MAG);
+        }
         this->p_chgmix->mix_reset();
+    }
 
     // mohan update 2012-06-05
     this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband();
diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
index 0bdf75045f..07779f62ec 100644
--- a/source/module_esolver/esolver_ks_pw.cpp
+++ b/source/module_esolver/esolver_ks_pw.cpp
@@ -492,9 +492,19 @@ void ESolver_KS_PW<T, Device>::othercalculation(const int istep)
 template <typename T, typename Device>
 void ESolver_KS_PW<T, Device>::eachiterinit(const int istep, const int iter)
 {
-    if (iter == 1)
+    if (iter == 1 || iter == GlobalV::MIXING_RESTART)
+    {
+        if (iter == GlobalV::MIXING_RESTART) // delete mixing and re-construct it to restart 
+        {
+            this->p_chgmix->set_mixing(GlobalV::MIXING_MODE,
+                                GlobalV::MIXING_BETA,
+                                GlobalV::MIXING_NDIM,
+                                GlobalV::MIXING_GG0,
+                                GlobalV::MIXING_TAU,
+                                GlobalV::MIXING_BETA_MAG);
+        }
         this->p_chgmix->mix_reset();
-
+    }
     // mohan move harris functional to here, 2012-06-05
     // use 'rho(in)' and 'v_h and v_xc'(in)
     this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband();
diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp
index 7a8d00a838..2f6b304a1c 100644
--- a/source/module_io/input.cpp
+++ b/source/module_io/input.cpp
@@ -303,6 +303,7 @@ void Input::Default(void)
     mixing_mode = "broyden";
     mixing_beta = -10;
     mixing_ndim = 8;
+    mixing_restart = 0;
     mixing_gg0 = 1.00;       // use Kerker defaultly
     mixing_beta_mag = -10.0; // only set when nspin == 2 || nspin == 4
     mixing_gg0_mag = 0.0;    // defaultly exclude Kerker from mixing magnetic density
@@ -1258,6 +1259,10 @@ bool Input::Read(const std::string& fn)
         {
             read_value(ifs, mixing_ndim);
         }
+        else if (strcmp("mixing_restart", word) == 0)
+        {
+            read_value(ifs, mixing_restart);
+        }
         else if (strcmp("mixing_gg0", word) == 0) // mohan add 2014-09-27
         {
             read_value(ifs, mixing_gg0);
@@ -3325,6 +3330,7 @@ void Input::Bcast()
     Parallel_Common::bcast_string(mixing_mode);
     Parallel_Common::bcast_double(mixing_beta);
     Parallel_Common::bcast_int(mixing_ndim);
+    Parallel_Common::bcast_int(mixing_restart);
     Parallel_Common::bcast_double(mixing_gg0); // mohan add 2014-09-27
     Parallel_Common::bcast_double(mixing_beta_mag);
     Parallel_Common::bcast_double(mixing_gg0_mag);
diff --git a/source/module_io/input.h b/source/module_io/input.h
index d0344b8944..ce723bab30 100644
--- a/source/module_io/input.h
+++ b/source/module_io/input.h
@@ -232,6 +232,7 @@ class Input
     std::string mixing_mode; // "plain","broyden",...
     double mixing_beta; // 0 : no_mixing
     int mixing_ndim; // used in Broyden method
+    int mixing_restart;
     double mixing_gg0; // used in kerker method. mohan add 2014-09-27
     double mixing_beta_mag;
     double mixing_gg0_mag;
diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp
index a52245d05c..8e00b1f774 100644
--- a/source/module_io/input_conv.cpp
+++ b/source/module_io/input_conv.cpp
@@ -750,6 +750,7 @@ void Input_Conv::Convert(void)
     GlobalV::MIXING_MODE = INPUT.mixing_mode;
     GlobalV::MIXING_BETA = INPUT.mixing_beta;
     GlobalV::MIXING_NDIM = INPUT.mixing_ndim;
+    GlobalV::MIXING_RESTART = INPUT.mixing_restart;
     GlobalV::MIXING_GG0 = INPUT.mixing_gg0;
     GlobalV::MIXING_BETA_MAG = INPUT.mixing_beta_mag;
     GlobalV::MIXING_GG0_MAG = INPUT.mixing_gg0_mag;
diff --git a/source/module_io/parameter_pool.cpp b/source/module_io/parameter_pool.cpp
index cf3e611ab3..524df9de87 100644
--- a/source/module_io/parameter_pool.cpp
+++ b/source/module_io/parameter_pool.cpp
@@ -818,6 +818,10 @@ void input_parameters_set(std::map<std::string, InputParameter> input_parameters
     {
         INPUT.mixing_ndim = *static_cast<int*>(input_parameters["mixing_ndim"].get());
     }
+    else if (input_parameters.count("mixing_restart") != 0)
+    {
+        INPUT.mixing_restart = *static_cast<int*>(input_parameters["mixing_restart"].get());
+    }
     else if (input_parameters.count("mixing_gg0") != 0)
     {
         INPUT.mixing_gg0 = *static_cast<double*>(input_parameters["mixing_gg0"].get());
diff --git a/source/module_io/test/input_conv_test.cpp b/source/module_io/test/input_conv_test.cpp
index f0d7e43f68..a566827792 100644
--- a/source/module_io/test/input_conv_test.cpp
+++ b/source/module_io/test/input_conv_test.cpp
@@ -183,6 +183,7 @@ TEST_F(InputConvTest, Conv)
 	EXPECT_EQ(GlobalV::sc_mag_switch,0);
     EXPECT_TRUE(GlobalV::decay_grad_switch);
     EXPECT_EQ(GlobalV::sc_file, "sc.json");
+	EXPECT_EQ(GlobalV::MIXING_RESTART,0);
 }
 
 TEST_F(InputConvTest, ConvRelax)
diff --git a/source/module_io/test/input_test_para.cpp b/source/module_io/test/input_test_para.cpp
index 7cc4d6bc73..d005fdfccc 100644
--- a/source/module_io/test/input_test_para.cpp
+++ b/source/module_io/test/input_test_para.cpp
@@ -381,6 +381,7 @@ TEST_F(InputParaTest, Bcast)
     EXPECT_TRUE(INPUT.mdp.dump_virial);
     EXPECT_FALSE(INPUT.mixing_tau);
     EXPECT_FALSE(INPUT.mixing_dftu);
+    EXPECT_EQ(INPUT.mixing_restart,0);
     EXPECT_EQ(INPUT.out_bandgap, 0);
     EXPECT_EQ(INPUT.out_mat_t, 0);
 
diff --git a/source/module_io/test/write_input_test.cpp b/source/module_io/test/write_input_test.cpp
index d61133715d..8dccb5627a 100644
--- a/source/module_io/test/write_input_test.cpp
+++ b/source/module_io/test/write_input_test.cpp
@@ -384,13 +384,16 @@ TEST_F(write_input, Mixing7)
     std::string output((std::istreambuf_iterator<char>(ifs)), std::istreambuf_iterator<char>());
     EXPECT_THAT(output, testing::HasSubstr("#Parameters (7.Charge Mixing)"));
     EXPECT_THAT(output, testing::HasSubstr("mixing_type                    broyden #plain; pulay; broyden"));
-    EXPECT_THAT(output,
-                testing::HasSubstr("mixing_beta                    0.7 #mixing parameter: 0 means no new charge"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_beta                    0.7 #mixing parameter: 0 means no new charge"));
     EXPECT_THAT(output, testing::HasSubstr("mixing_ndim                    8 #mixing dimension in pulay or broyden"));
     EXPECT_THAT(output, testing::HasSubstr("mixing_gg0                     0 #mixing parameter in kerker"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_beta_mag                -10 #mixing parameter for magnetic density"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_gg0_mag                 0 #mixing parameter in kerker"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_gg0_min                 0.1 #the minimum kerker coefficient"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_angle                   -10 #angle mixing parameter for non-colinear calculations"));
     EXPECT_THAT(output, testing::HasSubstr("mixing_tau                     0 #whether to mix tau in mGGA calculation"));
-    EXPECT_THAT(output,
-                testing::HasSubstr("mixing_dftu                    0 #whether to mix locale in DFT+U calculation"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_dftu                    0 #whether to mix locale in DFT+U calculation"));
+    EXPECT_THAT(output, testing::HasSubstr("mixing_restart                 0 #which step to restart mixing during SCF"));
     EXPECT_THAT(output, testing::HasSubstr(""));
     ifs.close();
     remove("write_input_test.log");
diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp
index 256a7dce73..e4c2c68464 100644
--- a/source/module_io/write_input.cpp
+++ b/source/module_io/write_input.cpp
@@ -248,6 +248,7 @@ ModuleBase::GlobalFunc::OUTP(ofs, "out_bandgap", out_bandgap, "if true, print ou
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_type", mixing_mode, "plain; pulay; broyden");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_beta", mixing_beta, "mixing parameter: 0 means no new charge");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_ndim", mixing_ndim, "mixing dimension in pulay or broyden");
+    ModuleBase::GlobalFunc::OUTP(ofs, "mixing_restart", mixing_restart, "which step to restart mixing during SCF");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_gg0", mixing_gg0, "mixing parameter in kerker");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_beta_mag", mixing_beta_mag, "mixing parameter for magnetic density");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_gg0_mag", mixing_gg0_mag, "mixing parameter in kerker");

From 3546506187e0ab5cbf7b6f30ca06fd79bc7f7736 Mon Sep 17 00:00:00 2001
From: jingan-181 <78459531+jingan-181@users.noreply.github.com>
Date: Thu, 25 Jan 2024 11:16:47 +0800
Subject: [PATCH 13/17] Update input-main.md (#3551)

Solve the format problem mentioned in issue 3543
---
 docs/advanced/input_files/input-main.md | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
index 42a4b3638a..8c675cab92 100644
--- a/docs/advanced/input_files/input-main.md
+++ b/docs/advanced/input_files/input-main.md
@@ -72,6 +72,7 @@
     - [mixing\_beta](#mixing_beta)
     - [mixing\_beta\_mag](#mixing_beta_mag)
     - [mixing\_ndim](#mixing_ndim)
+    - [mixing\_restart](#mixing_restart)
     - [mixing\_gg0](#mixing_gg0)
     - [mixing\_gg0\_mag](#mixing_gg0_mag)
     - [mixing\_gg0\_min](#mixing_gg0_min)
@@ -145,9 +146,9 @@
     - [out\_app\_flag](#out_app_flag)
     - [out\_ndigits](#out_ndigits)
     - [out\_interval](#out_interval)
-    - [out\_element\_info](#out_element_info)
     - [band\_print\_num](#band_print_num)
     - [bands\_to\_print](#bands_to_print)
+    - [out\_element\_info](#out_element_info)
     - [restart\_save](#restart_save)
     - [restart\_load](#restart_load)
     - [rpa](#rpa)
@@ -2801,9 +2802,9 @@ These variables are used to control berry phase and wannier90 interface paramete
 
 - **Type**: String
 - **Description**: the spin direction for the Wannier function calculation when nspin is set to 2
-  - "up": Calculate spin up for the Wannier function.
-  - "down": Calculate spin down for the Wannier function.
-- **Default**: "up"
+  - `up`: Calculate spin up for the Wannier function.
+  - `down`: Calculate spin down for the Wannier function.
+- **Default**: `up`
 
 ### out_wannier_mmn
 
@@ -2843,6 +2844,7 @@ These variables are used to control berry phase and wannier90 interface paramete
 - **Description**: write the "UNK.*" file in ASCII format or binary format.
   - 0: write the "UNK.*" file in binary format.
   - 1: write the "UNK.*" file in ASCII format (text file format).
+- **Default**: 1
 
 [back to top](#full-list-of-input-keywords)
 

From 2b254b5c0b56a09d4be32135bce4870e421f89d4 Mon Sep 17 00:00:00 2001
From: YI Zeping <18586016708@163.com>
Date: Thu, 25 Jan 2024 08:07:41 +0000
Subject: [PATCH 14/17] Build: fix compatibility issue against toolchain
 install (#3540)

* Fix for finding LibXC and ELPA

* For compatibility to previous routines

* syntax fix for FindELPA.cmake

* Update cmake/FindELPA.cmake

Co-authored-by: Chun Cai <amoycaic@gmail.com>

* Using CMake interface as default for finding LibXC

* update docs

* fix for FindLibxc: changing imcompatible if statement

* fix for FindLibxc: changing imcompatible if statement

* fix for FindLibxc: changing imcompatible if statement

* update docs for installing pkg-config

* Update FindLibxc.cmake

* Update FindLibxc.cmake

* remove previous LibXC routine in CMakeLists.txt

Co-authored-by: Chun Cai <amoycaic@gmail.com>

* Update easy_install.md with Makefile-built LibXC supported

* Update easy_install.md to include different behavior in different version on finding ELPA

* fix compatibility issue against toolchain

* Change default ELPA install routine to old one

---------

Co-authored-by: Chun Cai <amoycaic@gmail.com>
---
 cmake/FindELPA.cmake | 47 ++++++++++++++++++++++----------------------
 1 file changed, 24 insertions(+), 23 deletions(-)

diff --git a/cmake/FindELPA.cmake b/cmake/FindELPA.cmake
index 328b75035f..4105e47592 100644
--- a/cmake/FindELPA.cmake
+++ b/cmake/FindELPA.cmake
@@ -9,7 +9,26 @@
 
 find_package(PkgConfig)
 
-if(PKG_CONFIG_FOUND)
+find_path(ELPA_INCLUDE_DIRS
+    elpa/elpa.h
+    HINTS ${ELPA_DIR}
+    PATH_SUFFIXES "include" "include/elpa"
+    )
+if(USE_OPENMP)
+    find_library(ELPA_LINK_LIBRARIES
+    NAMES elpa_openmp elpa
+    HINTS ${ELPA_DIR}
+    PATH_SUFFIXES "lib"
+    )
+else()
+    find_library(ELPA_LINK_LIBRARIES
+    NAMES elpa
+    HINTS ${ELPA_DIR}
+    PATH_SUFFIXES "lib"
+    )
+endif()
+
+if(NOT ELPA_INCLUDE_DIRS AND PKG_CONFIG_FOUND)
   if(DEFINED ELPA_DIR)
     string(APPEND CMAKE_PREFIX_PATH ";${ELPA_DIR}")
   endif()
@@ -18,27 +37,9 @@ if(PKG_CONFIG_FOUND)
   else()
     pkg_search_module(ELPA REQUIRED IMPORTED_TARGET GLOBAL elpa)
   endif()
-else()
-  find_path(ELPA_INCLUDE_DIRS
-    elpa/elpa.h
-    HINTS ${ELPA_DIR}
-    PATH_SUFFIXES "include" "include/elpa"
-    )
-  if(USE_OPENMP)
-    find_library(ELPA_LINK_LIBRARIES
-      NAMES elpa_openmp elpa
-      HINTS ${ELPA_DIR}
-      PATH_SUFFIXES "lib"
-      )
-  else()
-    find_library(ELPA_LINK_LIBRARIES
-      NAMES elpa
-      HINTS ${ELPA_DIR}
-      PATH_SUFFIXES "lib"
-      )
-  endif()
-  #message(
-  #  "ELPA : We need pkg-config to get all information about the elpa library")
+elseif(NOT PKG_CONFIG_FOUND)
+  message(
+    "ELPA : We need pkg-config to get all information about the elpa library")
 endif()
 
 # Handle the QUIET and REQUIRED arguments and
@@ -48,7 +49,7 @@ find_package_handle_standard_args(ELPA DEFAULT_MSG ELPA_LINK_LIBRARIES ELPA_INCL
 
 # Copy the results to the output variables and target.
 if(ELPA_FOUND)
-    set(ELPA_LIBRARY ${ELPA_LINK_LIBRARIES})
+    list(GET ELPA_LINK_LIBRARIES 0 ELPA_LIBRARY)
     set(ELPA_INCLUDE_DIR ${ELPA_INCLUDE_DIRS})
 
     if(NOT TARGET ELPA::ELPA)

From 8dce30ac801bd6adb25f3be73b77ce347458217b Mon Sep 17 00:00:00 2001
From: Jie Li <76780849+jieli-matrix@users.noreply.github.com>
Date: Fri, 26 Jan 2024 12:01:09 +0800
Subject: [PATCH 15/17] Test: Configure performance tests for math libraries
 (#3511)

* add performace test of sphbes functions.

* fix benchmark cmake errors

* add dependencies for docker

* update docs

* add performance tests for sphbes

* add google benchmark

* rewrite benchmark tests in fixtures

* disable internal testing in benchmark

* merge benchmark into integration test

---------

Co-authored-by: StarGrys <771582678@qq.com>
---
 .github/workflows/test.yml                   |  2 +-
 CMakeLists.txt                               | 31 ++++++++-
 Dockerfile.cuda                              |  2 +-
 Dockerfile.gnu                               |  2 +-
 Dockerfile.intel                             |  2 +-
 docs/advanced/install.md                     | 10 +++
 docs/quick_start/easy_install.md             |  1 +
 source/module_base/test/CMakeLists.txt       |  7 ++
 source/module_base/test/perf_sphbes_test.cpp | 72 ++++++++++++++++++++
 9 files changed, 123 insertions(+), 6 deletions(-)
 create mode 100644 source/module_base/test/perf_sphbes_test.cpp

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 69d08d9b36..76f48347a8 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -29,7 +29,7 @@ jobs:
 
       - name: Build
         run: |
-          cmake -B build -DBUILD_TESTING=ON -DENABLE_DEEPKS=ON -DENABLE_LIBXC=ON -DENABLE_LIBRI=ON -DENABLE_PAW=ON
+          cmake -B build -DBUILD_TESTING=ON -DENABLE_DEEPKS=ON -DENABLE_LIBXC=ON -DENABLE_LIBRI=ON -DENABLE_PAW=ON -DENABLE_GOOGLEBENCH=ON
           cmake --build build -j8
           cmake --install build
 
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0ee1e3bbc0..e390578a2d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -34,6 +34,7 @@ option(DEBUG_INFO "Print message for developers to debug." OFF)
 option(ENABLE_NATIVE_OPTIMIZATION "Enable compilation optimization for the native machine's CPU type" OFF)
 option(COMMIT_INFO "Print commit information in log" ON)
 option(ENABLE_FFT_TWO_CENTER "Enable FFT-based two-center integral method." ON)
+option(ENABLE_GOOGLEBENCH "Enable GOOGLE-benchmark usage." OFF)
 option(ENABLE_RAPIDJSON "Enable rapid-json usage." OFF)
 
 
@@ -591,6 +592,25 @@ if(INFO)
   # modifications on blas_connector and lapack_connector
 endif()
 
+#  Add performance test in abacus
+IF (ENABLE_GOOGLEBENCH)
+  set(BUILD_TESTING ON)
+  find_package(benchmark HINTS ${BENCHMARK_DIR})
+  if(NOT ${benchmark_FOUND})
+    set(BENCHMARK_USE_BUNDLED_GTEST OFF)
+    include(FetchContent)
+    FetchContent_Declare(
+      benchmark
+      GIT_REPOSITORY https://github.com/google/benchmark.git
+      GIT_TAG "origin/main"
+      GIT_SHALLOW TRUE
+      GIT_PROGRESS TRUE 
+    )
+    set(BENCHMARK_ENABLE_TESTING OFF)
+    FetchContent_MakeAvailable(benchmark)
+  endif()
+endif()
+
 IF (BUILD_TESTING)
   set_if_higher(CMAKE_CXX_STANDARD 14) # Required in orbital
   include(CTest)
@@ -620,8 +640,14 @@ IF (BUILD_TESTING)
     endif()
 
     #dependencies & link library
-    target_link_libraries(${UT_TARGET} ${UT_LIBS}
-      Threads::Threads GTest::gtest_main GTest::gmock_main)
+    if(ENABLE_GOOGLEBENCH)
+      target_link_libraries(${UT_TARGET} ${UT_LIBS}
+        Threads::Threads GTest::gtest_main GTest::gmock_main benchmark::benchmark)
+      else()
+        target_link_libraries(${UT_TARGET} ${UT_LIBS}
+          Threads::Threads GTest::gtest_main GTest::gmock_main)
+    endif()
+      
     if(USE_OPENMP)
       target_link_libraries(${UT_TARGET} OpenMP::OpenMP_CXX)
     endif()
@@ -631,6 +657,7 @@ IF (BUILD_TESTING)
       WORKING_DIRECTORY $<TARGET_FILE_DIR:${UT_TARGET}>
     )
   endfunction(AddTest)
+
 endif()
 
 add_subdirectory(source)
diff --git a/Dockerfile.cuda b/Dockerfile.cuda
index 719f7c4278..e950f097f9 100644
--- a/Dockerfile.cuda
+++ b/Dockerfile.cuda
@@ -2,7 +2,7 @@ FROM nvidia/cuda:12.2.0-devel-ubuntu22.04
 
 RUN apt update && apt install -y --no-install-recommends \
     libopenblas-openmp-dev liblapack-dev libscalapack-mpi-dev libelpa-dev libfftw3-dev libcereal-dev \
-    libxc-dev libgtest-dev libgmock-dev python3-numpy \
+    libxc-dev libgtest-dev libgmock-dev libbenchmark-dev python3-numpy \
     bc cmake git g++ make bc time sudo unzip vim wget
 
 ENV GIT_SSL_NO_VERIFY=true TERM=xterm-256color \
diff --git a/Dockerfile.gnu b/Dockerfile.gnu
index 0b6b45d248..060d930563 100644
--- a/Dockerfile.gnu
+++ b/Dockerfile.gnu
@@ -1,7 +1,7 @@
 FROM ubuntu:22.04
 RUN apt update && apt install -y --no-install-recommends \
     libopenblas-openmp-dev liblapack-dev libscalapack-mpi-dev libelpa-dev libfftw3-dev libcereal-dev \
-    libxc-dev libgtest-dev libgmock-dev python3-numpy \
+    libxc-dev libgtest-dev libgmock-dev libbenchmark-dev python3-numpy \
     bc cmake git g++ make bc time sudo unzip vim wget gfortran
 
 ENV GIT_SSL_NO_VERIFY=true TERM=xterm-256color \
diff --git a/Dockerfile.intel b/Dockerfile.intel
index 6cac8c9f5f..3947f05b9e 100644
--- a/Dockerfile.intel
+++ b/Dockerfile.intel
@@ -2,7 +2,7 @@ FROM ubuntu:22.04
 
 RUN apt-get update && apt-get install -y \
     bc cmake git gnupg gcc g++ python3-numpy sudo wget vim unzip \
-    libcereal-dev libxc-dev libgtest-dev libgmock-dev
+    libcereal-dev libxc-dev libgtest-dev libgmock-dev libbenchmark-dev
 
 # Following steps by https://software.intel.com/content/www/us/en/develop/documentation/installation-guide-for-intel-oneapi-toolkits-linux/top/installation/install-using-package-managers/apt.html .
 RUN wget -O- https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB \
diff --git a/docs/advanced/install.md b/docs/advanced/install.md
index e929fac34c..d6201a060f 100644
--- a/docs/advanced/install.md
+++ b/docs/advanced/install.md
@@ -69,6 +69,16 @@ After building and installing, unit tests can be performed with `ctest`.
 
 To run a subset of unit test, use `ctest -R <test-match-pattern>` to perform tests with name matched by given pattern.
 
+## Build Performance Tests
+
+To build performance tests for ABACUS, define `ENABLE_GOOGLEBENCH` flag. You can also specify the path to a local installation of [Google Benchmark](https://github.com/google/benchmark.git) by setting `BENCHMARK_DIR` flags. If not found locally, the configuration process will try to download it automatically.
+
+```bash
+cmake -B build -DENABLE_GOOGLEBENCH=1
+```
+
+Google Benchmark requires Google Test to build and run the tests. When setting `ENABLE_GOOGLEBENCH` to ON, `BUILD_TESTING` is automatically enabled. After building and installing, performance tests can be executed with `ctest`.
+
 ## Build with CUDA support
 
 ### Extra prerequisites
diff --git a/docs/quick_start/easy_install.md b/docs/quick_start/easy_install.md
index f2173d569e..4089e303a3 100644
--- a/docs/quick_start/easy_install.md
+++ b/docs/quick_start/easy_install.md
@@ -126,6 +126,7 @@ Here, 'build' is the path for building ABACUS; and '-D' is used for setting up s
   - `ENABLE_LIBRI=OFF`: [Enable LibRI](../advanced/install.md#add-libri-support) to suppport variety of functionals. If `LIBRI_DIR` and `LIBCOMM_DIR` is defined, `ENABLE_LIBRI` will set to 'ON'.
   - `USE_OPENMP=ON`: Enable OpenMP support. Building ABACUS without OpenMP is not fully tested yet.
   - `BUILD_TESTING=OFF`: [Build unit tests](../advanced/install.md#build-unit-tests).
+  - `ENABLE_GOOGLEBENCH=OFF`: [Build performance tests](../advanced/install.md#build-performance-tests)
   - `ENABLE_MPI=ON`: Enable MPI parallel compilation. If set to `OFF`, a serial version of ABACUS with PW basis only will be compiled. Currently serial version of ABACUS with LCAO basis is not supported yet, so `ENABLE_LCAO` will be automatically set to `OFF`.
   - `ENABLE_COVERAGE=OFF`: Build ABACUS executable supporting [coverage analysis](../CONTRIBUTING.md#generating-code-coverage-report). This feature has a drastic impact on performance.
   - `ENABLE_ASAN=OFF`: Build with Address Sanitizer. This feature would help detecting memory problems.
diff --git a/source/module_base/test/CMakeLists.txt b/source/module_base/test/CMakeLists.txt
index 185727ef93..008df422e5 100644
--- a/source/module_base/test/CMakeLists.txt
+++ b/source/module_base/test/CMakeLists.txt
@@ -217,6 +217,13 @@ AddTest(
   SOURCES assoc_laguerre_test.cpp ../assoc_laguerre.cpp ../tool_quit.cpp ../global_variable.cpp ../global_file.cpp ../global_function.cpp ../memory.cpp ../timer.cpp
   LIBS ${math_libs} formatter
 )
+if(ENABLE_GOOGLEBENCH)
+  AddTest(
+    TARGET perf_sphbes
+    LIBS formatter
+    SOURCES perf_sphbes_test.cpp ../math_sphbes.cpp ../timer.cpp 
+  )
+endif()
 
 if(ENABLE_RAPIDJSON)
   AddTest(
diff --git a/source/module_base/test/perf_sphbes_test.cpp b/source/module_base/test/perf_sphbes_test.cpp
new file mode 100644
index 0000000000..4c574baa8e
--- /dev/null
+++ b/source/module_base/test/perf_sphbes_test.cpp
@@ -0,0 +1,72 @@
+#include"../math_sphbes.h"
+#include<fstream>
+#include <benchmark/benchmark.h>
+#include <iostream>
+#include <cstring>
+#include <cmath>
+
+/************************************************
+*  performace test of class Sphbes
+***********************************************/
+
+/**
+ * Tested function: 
+ *      - sphbesj
+ *      - Spherical_Bessel
+ */
+
+class PerfSphbes : public benchmark::Fixture {
+public:
+    const double q = 1;
+    const int n = 1000;
+    double stop = 1000.0;
+    double dr = 0.0;
+    double* rc, *rinf, *jc, *jinf;
+    void SetUp(const benchmark::State& state){
+        const double rcut = state.range(0) + 0.5;
+        rc = new double[n + 10]; 
+        rinf = new double[n + 10];
+        jc = new double[n + 10];
+        jinf = new double[n + 10];
+
+        // generate data points in (0, rcut] in log scale
+        double rmin = 0.0001;
+        double log_rmin = std::log(rmin);
+        double log_rcut = std::log(rcut);
+        dr = (log_rcut - log_rmin) / (n-1);
+        memset(rc, 0, (n+10) * sizeof(double));
+        for (int i = 0; i < n; i++)
+            rc[i] = std::exp(log_rmin + i * dr);
+        
+        // generate data points in [rcut, stop] in linear scale
+        memset(rinf, 0, (n+10) * sizeof(double));
+        rinf[0] = rcut;
+        dr = (stop - rcut) / (n-1);
+        for (int i = 1; i < n; i++)
+            rinf[i] += rinf[i-1] + dr;
+    }
+    void TearDown(const benchmark::State& state){
+        delete[] rc;
+        delete[] rinf;
+        delete[] jc;
+        delete[]  jinf;
+    }
+};    
+
+BENCHMARK_DEFINE_F(PerfSphbes, BM_Spherical_Bessel)(benchmark::State& state) {
+    for (auto _ : state) {
+        ModuleBase::Sphbes::Spherical_Bessel(n, rc, q, state.range(0), jc);
+        ModuleBase::Sphbes::Spherical_Bessel(n, rinf, q, state.range(0), jinf);
+    }
+}
+
+BENCHMARK_DEFINE_F(PerfSphbes, BM_sphbesj)(benchmark::State& state) {
+    for (auto _ : state) {
+        ModuleBase::Sphbes::sphbesj(n, rc, q, state.range(0), jc);
+        ModuleBase::Sphbes::sphbesj(n, rinf, q, state.range(0), jinf);
+    }
+}
+
+BENCHMARK_REGISTER_F(PerfSphbes, BM_sphbesj)->DenseRange(0, 11, 1)->Unit(benchmark::kMicrosecond);
+BENCHMARK_REGISTER_F(PerfSphbes, BM_Spherical_Bessel)->DenseRange(0, 11, 1)->Unit(benchmark::kMicrosecond);
+BENCHMARK_MAIN(); 
\ No newline at end of file

From 6a8ac74581b76a9750a065b38d1374cff259bdd0 Mon Sep 17 00:00:00 2001
From: Haozhi Han <haozhi.han@stu.pku.edu.cn>
Date: Sat, 27 Jan 2024 10:36:19 +0800
Subject: [PATCH 16/17] Update hsolver_pw.cpp (#3556)

when use_uspp==false, overlap matrix should be E.
---
 source/module_hsolver/hsolver_pw.cpp | 22 ++++++++++++++++++----
 1 file changed, 18 insertions(+), 4 deletions(-)

diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp
index 16fa5f335b..ae784d2009 100644
--- a/source/module_hsolver/hsolver_pw.cpp
+++ b/source/module_hsolver/hsolver_pw.cpp
@@ -624,17 +624,31 @@ void HSolverPW<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T, Device>* hm, psi::P
         hm->ops->hPsi(info);
         ModuleBase::timer::tick("DiagoCG_New", "hpsi_func");
     };
-    auto spsi_func = [hm](const ct::Tensor& psi_in, ct::Tensor& spsi_out) {
+    auto spsi_func = [this, hm](const ct::Tensor& psi_in, ct::Tensor& spsi_out) {
         ModuleBase::timer::tick("DiagoCG_New", "spsi_func");
         // psi_in should be a 2D tensor: 
         // psi_in.shape() = [nbands, nbasis]
         const auto ndim = psi_in.shape().ndim();
         REQUIRES_OK(ndim <= 2, "dims of psi_in should be less than or equal to 2");
-        // Convert a Tensor object to a psi::Psi object
-        hm->sPsi(psi_in.data<T>(), spsi_out.data<T>(), 
+
+        if (GlobalV::use_uspp)
+        {
+            // Convert a Tensor object to a psi::Psi object
+            hm->sPsi(psi_in.data<T>(), spsi_out.data<T>(), 
             ndim == 1 ? psi_in.NumElements() : psi_in.shape().dim_size(1), 
             ndim == 1 ? psi_in.NumElements() : psi_in.shape().dim_size(1), 
             ndim == 1 ? 1 : psi_in.shape().dim_size(0));
+        } else
+        {
+            psi::memory::synchronize_memory_op<T, Device, Device>()(
+                this->ctx,
+                this->ctx,
+                spsi_out.data<T>(),
+                psi_in.data<T>(),
+                static_cast<size_t>((ndim == 1 ? 1 : psi_in.shape().dim_size(0))
+                                    * (ndim == 1 ? psi_in.NumElements() : psi_in.shape().dim_size(1))));
+        }
+        
         ModuleBase::timer::tick("DiagoCG_New", "spsi_func");
     };
     auto psi_tensor = ct::TensorMap(
@@ -776,4 +790,4 @@ template class HSolverPW<std::complex<float>, psi::DEVICE_GPU>;
 template class HSolverPW<std::complex<double>, psi::DEVICE_GPU>;
 #endif
 
-} // namespace hsolver
\ No newline at end of file
+} // namespace hsolver

From 8239efb4845e5224252b79632daf5aaae871862f Mon Sep 17 00:00:00 2001
From: Chun Cai <amoycaic@gmail.com>
Date: Sat, 27 Jan 2024 10:42:35 +0800
Subject: [PATCH 17/17] Fix: cuda build target (#3276)

* Fix: cuda buid target

* Update CMakeLists.txt

---------

Co-authored-by: Denghui Lu <denghuilu@pku.edu.cn>
---
 CMakeLists.txt | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index e390578a2d..9acc91959d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -429,8 +429,7 @@ endif()
 
 if(ENABLE_DEEPKS)
   # Torch uses outdated components to detech CUDA arch, causing failure on latest CUDA kits.
-  # See above for setting CMAKE_CUDA_ARCHITECTURES
-  set(TORCH_CUDA_ARCH_LIST CMAKE_CUDA_ARCHITECTURES)
+  # Set CMake variable TORCH_CUDA_ARCH_LIST in the form of "major.minor" if required.
   find_package(Torch REQUIRED)
   if(NOT Torch_VERSION VERSION_LESS "2.1.0")
     set_if_higher(CMAKE_CXX_STANDARD 17)