diff --git a/python/pyabacus/src/ModuleNAO/py_m_nao.cpp b/python/pyabacus/src/ModuleNAO/py_m_nao.cpp
index 45e6df2710..22db0f99b0 100644
--- a/python/pyabacus/src/ModuleNAO/py_m_nao.cpp
+++ b/python/pyabacus/src/ModuleNAO/py_m_nao.cpp
@@ -82,9 +82,9 @@ void bind_m_nao(py::module& m)
         .def("symbol", &RadialCollection::symbol, "itype"_a)
         .def_property_readonly("ntype", &RadialCollection::ntype)
         .def("lmax", overload_cast_<const int>()(&RadialCollection::lmax, py::const_), "itype"_a)
-        .def_property_readonly("lmax", overload_cast_<>()(&RadialCollection::lmax, py::const_))
+        .def("lmax", overload_cast_<>()(&RadialCollection::lmax, py::const_))
         .def("rcut_max", overload_cast_<const int>()(&RadialCollection::rcut_max, py::const_), "itype"_a)
-        .def_property_readonly("rcut_max", overload_cast_<>()(&RadialCollection::rcut_max, py::const_))
+        .def("rcut_max", overload_cast_<>()(&RadialCollection::rcut_max, py::const_))
         .def("nzeta", &RadialCollection::nzeta, "itype"_a, "l"_a)
         .def("nzeta_max", overload_cast_<const int>()(&RadialCollection::nzeta_max, py::const_), "itype"_a)
         .def("nzeta_max", overload_cast_<>()(&RadialCollection::nzeta_max, py::const_))
diff --git a/python/pyabacus/src/pyabacus/ModuleNAO/_module_nao.py b/python/pyabacus/src/pyabacus/ModuleNAO/_module_nao.py
index b8ed03a1a4..132a955cbc 100644
--- a/python/pyabacus/src/pyabacus/ModuleNAO/_module_nao.py
+++ b/python/pyabacus/src/pyabacus/ModuleNAO/_module_nao.py
@@ -83,19 +83,21 @@ def symbol(self, itype: int) -> str:
     def ntype(self) -> int: 
         return super().ntype
     
-    def lmax(self, itype: int) -> int: 
-        return super().lmax(itype)
-    
-    @property
-    def lmax(self) -> int: 
-        return super().lmax
-    
-    def rcut_max(self, itype: int) -> float: 
-        return super().rcut_max(itype)
+    @overload
+    def lmax(self) -> int: ...
+    @overload
+    def lmax(self, itype: int) -> int: ...
 
-    @property
-    def rcut_max(self) -> float: 
-        return super().rcut_max
+    def lmax(self, *args, **kwargs): 
+        return super().lmax(*args, **kwargs)
+
+    @overload
+    def rcut_max(self) -> float: ...
+    @overload
+    def rcut_max(self, itype: int) -> float: ...
+
+    def rcut_max(self, *args, **kwargs): 
+        return super().rcut_max(*args, **kwargs)
     
     def nzeta(self, itype: int, l: int) -> int: 
         return super().nzeta(itype, l)
diff --git a/python/pyabacus/src/pyabacus/io/stru.py b/python/pyabacus/src/pyabacus/io/stru.py
index 1b337e74a7..0327898fc4 100644
--- a/python/pyabacus/src/pyabacus/io/stru.py
+++ b/python/pyabacus/src/pyabacus/io/stru.py
@@ -107,6 +107,7 @@ def read_stru(fpath: str) -> Dict[str, Any]:
     """Read an ABACUS STRU file and return its content as a dictionary."""
     block_title = ['ATOMIC_SPECIES',
                    'NUMERICAL_ORBITAL',
+                   'NUMERICAL_DESCRIPTOR',
                    'LATTICE_CONSTANT',
                    'LATTICE_PARAMETER',
                    'LATTICE_VECTORS',
@@ -144,6 +145,10 @@ def _trim(line):
         for i, s in enumerate(stru['species']):
             s['orb_file'] = blocks['NUMERICAL_ORBITAL'][i].strip()
 
+    #============ NUMERICAL_DESCRIPTOR ============
+    if 'NUMERICAL_DESCRIPTOR' in blocks:
+        stru['desc'] = blocks['NUMERICAL_DESCRIPTOR'][0].strip()
+
     #============ ATOMIC_POSITIONS ============
     stru['coord_type'] = blocks['ATOMIC_POSITIONS'][0]
     index = {s['symbol']: i for i, s in enumerate(stru['species'])}
diff --git a/python/pyabacus/tests/test_m_nao.py b/python/pyabacus/tests/test_m_nao.py
index d103c7fe2c..9a13d1b3fd 100644
--- a/python/pyabacus/tests/test_m_nao.py
+++ b/python/pyabacus/tests/test_m_nao.py
@@ -50,8 +50,8 @@ def test_rc():
     assert orb.symbol(3) == 'Fe'
 
     assert orb.ntype == 4
-    assert orb.lmax == 3
-    assert orb.rcut_max == 10.0
+    assert orb.lmax() == 3
+    assert orb.rcut_max() == 10.0
 
     assert orb.nzeta(0,0) == 2
     assert orb.nzeta(1,0) == 2
@@ -85,7 +85,7 @@ def test_twocenterintegrator():
     alpha.build(1, [file_list[0]])
 
     dr = 0.01  # R spacing
-    rmax = max(orb.rcut_max, alpha.rcut_max)
+    rmax = max(orb.rcut_max(), alpha.rcut_max())
     cutoff = 2.0 * rmax
     nr = int(rmax / dr) + 1
 
diff --git a/source/source_basis/module_nao/two_center_integrator.cpp b/source/source_basis/module_nao/two_center_integrator.cpp
index ebc7f71141..295946bb9e 100644
--- a/source/source_basis/module_nao/two_center_integrator.cpp
+++ b/source/source_basis/module_nao/two_center_integrator.cpp
@@ -48,6 +48,12 @@ void TwoCenterIntegrator::calculate(const int itype1,
         return;
     }
 
+    if (m1 > l1 || m1 < -l1 || m2 > l2 || m2 < -l2)
+    {
+        ModuleBase::WARNING("TwoCenterIntegrator", "m should be in range [-l, l].");
+        return;
+    }
+
     // unit vector along R
     ModuleBase::Vector3<double> uR = (R == 0.0 ? ModuleBase::Vector3<double>(0., 0., 1.) : vR / R);