diff --git a/source/source_cell/module_neighlist/CMakeLists.txt b/source/source_cell/module_neighlist/CMakeLists.txt
index 6f54f71513f..b40a0d4de29 100644
--- a/source/source_cell/module_neighlist/CMakeLists.txt
+++ b/source/source_cell/module_neighlist/CMakeLists.txt
@@ -3,6 +3,8 @@ add_library(
     OBJECT
     bin_manager.cpp
     neighbor_search.cpp
+    page_allocator.cpp
+    unitcell_lite.cpp
 )
 
 if(ENABLE_COVERAGE)
diff --git a/source/source_cell/module_neighlist/atom_provider.h b/source/source_cell/module_neighlist/atom_provider.h
new file mode 100644
index 00000000000..3087148dac5
--- /dev/null
+++ b/source/source_cell/module_neighlist/atom_provider.h
@@ -0,0 +1,74 @@
+#ifndef ATOM_PROVIDER_H
+#define ATOM_PROVIDER_H
+
+#include "source_base/vector3.h"
+#include "source_base/matrix3.h"
+
+/**
+ * @brief Interface for providing atom and lattice information.
+ *
+ * This abstract interface defines the minimum set of methods needed by
+ * the neighbor search module to access atom positions and lattice parameters.
+ * Any class implementing this interface can be used with NeighborSearch.
+ *
+ * @see UnitCell
+ * @see UnitCellLite
+ */
+class AtomProvider
+{
+public:
+    /**
+     * @brief Default destructor.
+     */
+    virtual ~AtomProvider() = default;
+
+    /**
+     * @brief Get the lattice constant.
+     * @return Lattice constant in Bohr.
+     */
+    virtual double get_lat0() const = 0;
+
+    /**
+     * @brief Get the volume of the unit cell.
+     * @return Unit cell volume in Bohr^3.
+     */
+    virtual double get_omega() const = 0;
+
+    /**
+     * @brief Get the lattice vectors.
+     * @return Const reference to the 3x3 lattice vector matrix.
+     */
+    virtual const ModuleBase::Matrix3& get_latvec() const = 0;
+
+    /**
+     * @brief Get the total number of atoms.
+     * @return Total atom count.
+     */
+    virtual int get_natom() const = 0;
+
+    /**
+     * @brief Get the number of atoms of a specific type.
+     * @param i Type index.
+     * @return Number of atoms of type i.
+     */
+    virtual int get_na(int i) const = 0;
+
+    /**
+     * @brief Get the number of atom types.
+     * @return Number of atom types.
+     */
+    virtual int get_ntype() const = 0;
+
+    /**
+     * @brief Get the Cartesian coordinates of a specific atom.
+     *
+     * Returns the position of the j-th atom of type i.
+     *
+     * @param i Type index.
+     * @param j Atom index within type i.
+     * @return Cartesian position vector.
+     */
+    virtual ModuleBase::Vector3<double> get_tau(int i, int j) const = 0;
+};
+
+#endif // ATOM_PROVIDER_H
\ No newline at end of file
diff --git a/source/source_cell/module_neighlist/bin_manager.cpp b/source/source_cell/module_neighlist/bin_manager.cpp
index 2e64b49a881..1077b91dd75 100644
--- a/source/source_cell/module_neighlist/bin_manager.cpp
+++ b/source/source_cell/module_neighlist/bin_manager.cpp
@@ -1,9 +1,67 @@
 #include <limits>
 #include <cmath>
 #include <algorithm>
+#include <cassert>
 #include "bin_manager.h"
 
+// ========== Bin class implementation ==========
 
+int Bin::get_id_x() const {
+    return id_x_;
+}
+
+int Bin::get_id_y() const {
+    return id_y_;
+}
+
+int Bin::get_id_z() const {
+    return id_z_;
+}
+
+const std::vector<NeighborAtom>& Bin::get_atoms() const {
+    return atoms_;
+}
+
+void Bin::set_id(int ix, int iy, int iz) {
+    id_x_ = ix;
+    id_y_ = iy;
+    id_z_ = iz;
+}
+
+void Bin::clear_atoms() {
+    atoms_.clear();
+}
+
+void Bin::add_atom(const NeighborAtom& atom) {
+    atoms_.push_back(atom);
+}
+
+// ========== BinManager getter methods ==========
+
+int BinManager::get_nbinx() const {
+    return nbinx_;
+}
+
+int BinManager::get_nbiny() const {
+    return nbiny_;
+}
+
+int BinManager::get_nbinz() const {
+    return nbinz_;
+}
+
+int BinManager::get_total_bins() const {
+    return static_cast<int>(bins_.size());
+}
+
+int BinManager::get_bin_atom_count(int bin_index) const {
+    if (bin_index < 0 || bin_index >= static_cast<int>(bins_.size())) {
+        return 0;
+    }
+    return static_cast<int>(bins_[bin_index].get_atoms().size());
+}
+
+// ========== BinManager main methods ==========
 
 void BinManager::init_bins(
     double sr,
@@ -11,69 +69,63 @@ void BinManager::init_bins(
     const std::vector<NeighborAtom>& ghost_atoms
 )
 {
-    sradius = sr;
+    sradius_ = sr;
     if(inside_atoms.empty() && ghost_atoms.empty())
     {
-        x_min=y_min=z_min=0;
-        x_max=y_max=z_max=0;
-        nbinx=nbiny=nbinz=1;
-        bins.clear();
-        bins.resize(1);
+        x_min_ = y_min_ = z_min_ = 0;
+        x_max_ = y_max_ = z_max_ = 0;
+        nbinx_ = nbiny_ = nbinz_ = 1;
+        bins_.clear();
+        bins_.resize(1);
         return;
     }
 
-    x_min = y_min = z_min = std::numeric_limits<double>::max();
-
-    x_max = y_max = z_max = std::numeric_limits<double>::lowest();
-    
+    x_min_ = y_min_ = z_min_ = std::numeric_limits<double>::max();
+    x_max_ = y_max_ = z_max_ = std::numeric_limits<double>::lowest();
 
     auto update_bounds = [&](const std::vector<NeighborAtom>& atoms)
     {
         for (const auto& atom : atoms)
         {
-            x_min = std::min(x_min, atom.position_x);
-            x_max = std::max(x_max, atom.position_x);
+            x_min_ = std::min(x_min_, atom.position_x);
+            x_max_ = std::max(x_max_, atom.position_x);
 
-            y_min = std::min(y_min, atom.position_y);
-            y_max = std::max(y_max, atom.position_y);
+            y_min_ = std::min(y_min_, atom.position_y);
+            y_max_ = std::max(y_max_, atom.position_y);
 
-            z_min = std::min(z_min, atom.position_z);
-            z_max = std::max(z_max, atom.position_z);
+            z_min_ = std::min(z_min_, atom.position_z);
+            z_max_ = std::max(z_max_, atom.position_z);
         }
     };
 
     update_bounds(inside_atoms);
     update_bounds(ghost_atoms);
 
-    bin_sizex = bin_sizey = bin_sizez = sradius;
+    bin_sizex_ = bin_sizey_ = bin_sizez_ = sradius_;
 
-    nbinx = std::ceil((x_max - x_min) / bin_sizex);
-    nbiny = std::ceil((y_max - y_min) / bin_sizey);
-    nbinz = std::ceil((z_max - z_min) / bin_sizez);
+    nbinx_ = std::ceil((x_max_ - x_min_) / bin_sizex_);
+    nbiny_ = std::ceil((y_max_ - y_min_) / bin_sizey_);
+    nbinz_ = std::ceil((z_max_ - z_min_) / bin_sizez_);
 
-    nbinx = std::max(1, nbinx);
-    nbiny = std::max(1, nbiny);
-    nbinz = std::max(1, nbinz);
+    nbinx_ = std::max(1, nbinx_);
+    nbiny_ = std::max(1, nbiny_);
+    nbinz_ = std::max(1, nbinz_);
 
-    int nbins = nbinx * nbiny * nbinz;
+    int nbins = nbinx_ * nbiny_ * nbinz_;
 
-    bins.clear();
+    bins_.clear();
+    bins_.resize(nbins);
 
-    bins.resize(nbins);
-
-    for (int ix = 0; ix < nbinx; ++ix)
+    for (int ix = 0; ix < nbinx_; ++ix)
     {
-        for (int iy = 0; iy < nbiny; ++iy)
+        for (int iy = 0; iy < nbiny_; ++iy)
         {
-            for (int iz = 0; iz < nbinz; ++iz)
+            for (int iz = 0; iz < nbinz_; ++iz)
             {
-                int idx = ix * nbiny * nbinz + iy * nbinz + iz;
-
-                bins[idx].id_x = ix;
-                bins[idx].id_y = iy;
-                bins[idx].id_z = iz;
+                int idx = bin_index(ix, iy, iz);
 
-                bins[idx].atoms.clear();
+                bins_[idx].set_id(ix, iy, iz);
+                bins_[idx].clear_atoms();
             }
         }
     }
@@ -87,58 +139,63 @@ void BinManager::do_binning(
     auto bin_atom = [&](const NeighborAtom& atom)
     {
         int ix = std::min(
-            std::max(int((atom.position_x - x_min) / bin_sizex), 0),
-            nbinx - 1
+            std::max(int((atom.position_x - x_min_) / bin_sizex_), 0),
+            nbinx_ - 1
         );
 
         int iy = std::min(
-            std::max(int((atom.position_y - y_min) / bin_sizey), 0),
-            nbiny - 1
+            std::max(int((atom.position_y - y_min_) / bin_sizey_), 0),
+            nbiny_ - 1
         );
 
         int iz = std::min(
-            std::max(int((atom.position_z - z_min) / bin_sizez), 0),
-            nbinz - 1
+            std::max(int((atom.position_z - z_min_) / bin_sizez_), 0),
+            nbinz_ - 1
         );
 
-        int idx = ix * nbiny * nbinz + iy * nbinz + iz;
+        int idx = bin_index(ix, iy, iz);
 
-        bins[idx].atoms.push_back(atom);
+        bins_[idx].add_atom(atom);
     };
 
     for (const auto& atom : inside_atoms) bin_atom(atom);
-
     for (const auto& atom : ghost_atoms) bin_atom(atom);
 }
 
+int BinManager::bin_index(int ix, int iy, int iz) const {
+    return ix * nbiny_ * nbinz_ + iy * nbinz_ + iz;
+}
+
 void BinManager::build_atom_neighbors(
     NeighborList& neighbor_list,
     std::vector<NeighborAtom>& atoms
 )
 {
-    assert(atoms.size() == neighbor_list.numneigh.size());
+    assert(atoms.size() == static_cast<size_t>(neighbor_list.get_nlocal()));
 
-    double sradius2 = sradius * sradius;
+    double sradius2 = sradius_ * sradius_;
 
     neighbor_list.reset();
 
+    std::vector<int> neigh_tmp;
+
     for (int i = 0; i < atoms.size(); i++)
     {
-        std::vector<int> neigh_tmp;
+        neigh_tmp.clear();
 
         int ix = std::min(
-            std::max(int((atoms[i].position_x - x_min) / bin_sizex), 0),
-            nbinx - 1
+            std::max(int((atoms[i].position_x - x_min_) / bin_sizex_), 0),
+            nbinx_ - 1
         );
 
         int iy = std::min(
-            std::max(int((atoms[i].position_y - y_min) / bin_sizey), 0),
-            nbiny - 1
+            std::max(int((atoms[i].position_y - y_min_) / bin_sizey_), 0),
+            nbiny_ - 1
         );
 
         int iz = std::min(
-            std::max(int((atoms[i].position_z - z_min) / bin_sizez), 0),
-            nbinz - 1
+            std::max(int((atoms[i].position_z - z_min_) / bin_sizez_), 0),
+            nbinz_ - 1
         );
 
         for (int dx = -1; dx <= 1; dx++)
@@ -151,14 +208,14 @@ void BinManager::build_atom_neighbors(
                     int jy = iy + dy;
                     int jz = iz + dz;
 
-                    if (jx < 0 || jx >= nbinx ||
-                        jy < 0 || jy >= nbiny ||
-                        jz < 0 || jz >= nbinz)
+                    if (jx < 0 || jx >= nbinx_ ||
+                        jy < 0 || jy >= nbiny_ ||
+                        jz < 0 || jz >= nbinz_)
                         continue;
 
-                    int nidx = jx * nbiny * nbinz + jy * nbinz + jz;
+                    int nidx = bin_index(jx, jy, jz);
 
-                    for (const NeighborAtom& natom : bins[nidx].atoms)
+                    for (const NeighborAtom& natom : bins_[nidx].get_atoms())
                     {
                         double dx = atoms[i].position_x - natom.position_x;
                         double dy = atoms[i].position_y - natom.position_y;
@@ -173,30 +230,29 @@ void BinManager::build_atom_neighbors(
                     }
                 }
             }
-        } 
+        }
+
         int n = neigh_tmp.size();
 
-        //std::cout<<n<<std::endl;
+        int* ptr = neighbor_list.allocator_.allocate(n);
 
-        int* ptr = neighbor_list.allocator.allocate(n);
-    
         for (int k = 0; k < n; k++)
         {
             assert(ptr != nullptr);
             ptr[k] = neigh_tmp[k];
         }
 
-        neighbor_list.firstneigh[i] = ptr;
-        neighbor_list.numneigh[i] = n;
+        neighbor_list.firstneigh_[i] = ptr;
+        neighbor_list.numneigh_[i] = n;
     }
 }
 
 void BinManager::clear()
 {
-    for (auto& bin : bins)
+    for (auto& bin : bins_)
     {
-        bin.atoms.clear();
+        bin.clear_atoms();
     }
 
-    bins.clear();
+    bins_.clear();
 }
\ No newline at end of file
diff --git a/source/source_cell/module_neighlist/bin_manager.h b/source/source_cell/module_neighlist/bin_manager.h
index e8d3f0e9546..22b94d394a0 100644
--- a/source/source_cell/module_neighlist/bin_manager.h
+++ b/source/source_cell/module_neighlist/bin_manager.h
@@ -5,62 +5,219 @@
 #include "source_cell/module_neighlist/neighbor_atom.h"
 #include "source_cell/module_neighlist/neighbor_list.h"
 
-
-
+/**
+ * @brief A single bin in the 3D binning grid for neighbor search.
+ *
+ * Each bin stores atoms that fall within its spatial region,
+ * along with its position indices in the 3D grid.
+ */
 class Bin
 {
 public:
-    Bin()=default;
-    ~Bin()=default;
-
-    int id_x;
-    int id_y;
-    int id_z;
-
-    std::vector<NeighborAtom> atoms;
+    /**
+     * @brief Default constructor.
+     */
+    Bin() = default;
+
+    /**
+     * @brief Default destructor.
+     */
+    ~Bin() = default;
+
+    // ========== Getter methods ==========
+
+    /**
+     * @brief Get the X index of this bin in the grid.
+     * @return X index.
+     */
+    int get_id_x() const;
+
+    /**
+     * @brief Get the Y index of this bin in the grid.
+     * @return Y index.
+     */
+    int get_id_y() const;
+
+    /**
+     * @brief Get the Z index of this bin in the grid.
+     * @return Z index.
+     */
+    int get_id_z() const;
+
+    /**
+     * @brief Get the atoms stored in this bin.
+     * @return Const reference to the atom vector.
+     */
+    const std::vector<NeighborAtom>& get_atoms() const;
+
+    // ========== Setter methods (internal use) ==========
+
+    /**
+     * @brief Set the grid indices for this bin.
+     * @param ix X index.
+     * @param iy Y index.
+     * @param iz Z index.
+     */
+    void set_id(int ix, int iy, int iz);
+
+    /**
+     * @brief Clear all atoms from this bin.
+     */
+    void clear_atoms();
+
+    /**
+     * @brief Add an atom to this bin.
+     * @param atom The atom to add.
+     */
+    void add_atom(const NeighborAtom& atom);
+
+private:
+    /// X index in the 3D bin grid
+    int id_x_ = 0;
+
+    /// Y index in the 3D bin grid
+    int id_y_ = 0;
+
+    /// Z index in the 3D bin grid
+    int id_z_ = 0;
+
+    /// Atoms contained in this bin
+    std::vector<NeighborAtom> atoms_;
 };
 
-
+/**
+ * @brief Manager for binning atoms to accelerate neighbor search.
+ *
+ * This class implements a spatial binning strategy where atoms are
+ * sorted into a 3D grid of bins. Neighbor search then only needs to
+ * check adjacent bins, significantly reducing the search complexity.
+ *
+ * The workflow is:
+ * 1. Call init_bins() to set up the bin grid based on atom positions
+ * 2. Call do_binning() to assign atoms to bins
+ * 3. Call build_atom_neighbors() to construct the neighbor list
+ * 4. Call clear() to reset for the next search
+ */
 class BinManager
 {
 public:
-
+    /**
+     * @brief Initialize the bin grid based on atom positions and search radius.
+     *
+     * Computes the spatial bounds from atom positions and divides the
+     * region into bins of size approximately equal to the search radius.
+     *
+     * @param sr Search radius in lattice units.
+     * @param inside_atoms Atoms inside the local MPI domain.
+     * @param ghost_atoms Ghost atoms from neighboring domains.
+     */
     void init_bins(
         double sr,
         const std::vector<NeighborAtom>& inside_atoms,
         const std::vector<NeighborAtom>& ghost_atoms
     );
 
-
+    /**
+     * @brief Assign atoms to their corresponding bins.
+     *
+     * Must be called after init_bins(). Each atom is placed into the
+     * bin that contains its spatial position.
+     *
+     * @param inside_atoms Atoms inside the local MPI domain.
+     * @param ghost_atoms Ghost atoms from neighboring domains.
+     */
     void do_binning(
         const std::vector<NeighborAtom>& inside_atoms,
         const std::vector<NeighborAtom>& ghost_atoms
     );
 
-
+    /**
+     * @brief Build neighbor list by searching adjacent bins.
+     *
+     * For each atom, searches its containing bin and all adjacent bins
+     * to find neighbors within the search radius.
+     *
+     * @param neighbor_list Output neighbor list to populate.
+     * @param atoms Atoms for which to build neighbors.
+     */
     void build_atom_neighbors(
         NeighborList& neighbor_list,
         std::vector<NeighborAtom>& atoms
     );
 
-
+    /**
+     * @brief Clear all bins and reset internal state.
+     */
     void clear();
 
-
-    double sradius;
-
-    double x_min, y_min, z_min;
-    double x_max, y_max, z_max;
-
-    double bin_sizex;
-    double bin_sizey;
-    double bin_sizez;
-
-    int nbinx;
-    int nbiny;
-    int nbinz;
-
-    std::vector<Bin> bins;
+    // ========== Getter methods ==========
+
+    /**
+     * @brief Get the number of bins in X direction.
+     * @return Number of bins in X.
+     */
+    int get_nbinx() const;
+
+    /**
+     * @brief Get the number of bins in Y direction.
+     * @return Number of bins in Y.
+     */
+    int get_nbiny() const;
+
+    /**
+     * @brief Get the number of bins in Z direction.
+     * @return Number of bins in Z.
+     */
+    int get_nbinz() const;
+
+    /**
+     * @brief Get the total number of bins.
+     * @return Total bin count (nbinx * nbiny * nbinz).
+     */
+    int get_total_bins() const;
+
+    /**
+     * @brief Get the number of atoms in a specific bin.
+     * @param bin_index Index of the bin in the flat array.
+     * @return Number of atoms in that bin.
+     */
+    int get_bin_atom_count(int bin_index) const;
+
+private:
+    /// Search radius in lattice units
+    double sradius_ = 0.0;
+
+    /// Minimum coordinates of the binned region
+    double x_min_ = 0.0;
+    double y_min_ = 0.0;
+    double z_min_ = 0.0;
+
+    /// Maximum coordinates of the binned region
+    double x_max_ = 0.0;
+    double y_max_ = 0.0;
+    double z_max_ = 0.0;
+
+    /// Size of each bin in each direction
+    double bin_sizex_ = 0.0;
+    double bin_sizey_ = 0.0;
+    double bin_sizez_ = 0.0;
+
+    /// Number of bins in each direction
+    int nbinx_ = 1;
+    int nbiny_ = 1;
+    int nbinz_ = 1;
+
+    /// All bins in the 3D grid (stored as flat array)
+    std::vector<Bin> bins_;
+
+    /**
+     * @brief Compute the flat index for a bin from its 3D coordinates.
+     * @param ix X index of the bin.
+     * @param iy Y index of the bin.
+     * @param iz Z index of the bin.
+     * @return Flat index in the bins_ array.
+     */
+    int bin_index(int ix, int iy, int iz) const;
 };
 
 #endif // BIN_MANAGER_H
\ No newline at end of file
diff --git a/source/source_cell/module_neighlist/neighbor_atom.h b/source/source_cell/module_neighlist/neighbor_atom.h
index ed8d99c161a..9e9525c9e97 100644
--- a/source/source_cell/module_neighlist/neighbor_atom.h
+++ b/source/source_cell/module_neighlist/neighbor_atom.h
@@ -3,30 +3,90 @@
 
 #include <vector>
 
+/**
+ * @brief Represents an atom with neighbor search related properties.
+ *
+ * This class stores atom position, type, and index information used during
+ * neighbor search operations. It distinguishes between atoms inside the
+ * local MPI domain and ghost atoms from neighboring domains.
+ */
 class NeighborAtom
 {
 public:
+    /// X coordinate of the atom in Cartesian coordinates
     double position_x;
+
+    /// Y coordinate of the atom in Cartesian coordinates
     double position_y;
+
+    /// Z coordinate of the atom in Cartesian coordinates
     double position_z;
+
+    /// Atom type index
     int atom_type;
+
+    /// Index of the atom within its type
     int atom_index;
+
+    /// Unique atom ID across all domains and periodic images
     int atom_id;
-    //bool isghost;
+
+    /// Whether this atom is inside the local MPI domain
     bool is_inside;
 
+    /**
+     * @brief Construct a NeighborAtom.
+     *
+     * @param x X coordinate.
+     * @param y Y coordinate.
+     * @param z Z coordinate.
+     * @param type Atom type index.
+     * @param index Index within the atom type.
+     * @param id Unique atom ID.
+     */
     NeighborAtom(double x, double y, double z, int type, int index, int id)
         : position_x(x), position_y(y), position_z(z),
-          atom_type(type), atom_index(index), atom_id(id) {}
+          atom_type(type), atom_index(index), atom_id(id), is_inside(false) {}
 };
 
+/**
+ * @brief Input structure for neighbor search initialization.
+ *
+ * Contains atom data and spatial bounds computed from input atoms,
+ * used to initialize the binning grid.
+ */
 class InputAtoms
 {
 public:
+    /// List of input atoms
     std::vector<NeighborAtom> InputAtom;
-    double x_low, x_high, y_low, y_high, z_low, z_high;
+
+    /// Minimum X coordinate of the atom bounding box
+    double x_low;
+
+    /// Maximum X coordinate of the atom bounding box
+    double x_high;
+
+    /// Minimum Y coordinate of the atom bounding box
+    double y_low;
+
+    /// Maximum Y coordinate of the atom bounding box
+    double y_high;
+
+    /// Minimum Z coordinate of the atom bounding box
+    double z_low;
+
+    /// Maximum Z coordinate of the atom bounding box
+    double z_high;
+
+    /// Total number of atoms
     int n_atoms;
 
+    /**
+     * @brief Default constructor.
+     *
+     * Initializes bounds to zero and atom count to zero.
+     */
     InputAtoms()
         : x_low(0), x_high(0), y_low(0), y_high(0), z_low(0), z_high(0), n_atoms(0) {}
 };
diff --git a/source/source_cell/module_neighlist/neighbor_list.h b/source/source_cell/module_neighlist/neighbor_list.h
index 9ed0a4dc409..c14c80535f0 100644
--- a/source/source_cell/module_neighlist/neighbor_list.h
+++ b/source/source_cell/module_neighlist/neighbor_list.h
@@ -1,100 +1,42 @@
-#pragma once
+#ifndef NEIGHBOR_LIST_H
+#define NEIGHBOR_LIST_H
 
 #include <vector>
-#include<iostream>
-#include <cmath>
-#include <cassert>
+#include "page_allocator.h"
 
-class PageAllocator {
+class NeighborList
+{
 public:
-    struct Page {
-        std::vector<int> data;
-        int capacity;
-        int offset;
-    };
-
-    std::vector<Page> pages;
-    int pgsize;
-
-    static constexpr int DEFAULT_PGSIZE = 1024;
-
-    // Default constructor
-    PageAllocator() : pgsize(DEFAULT_PGSIZE) {
-        if (pgsize > 0) new_page();
-    }
+    NeighborList() = default;
+    ~NeighborList() = default;
 
-    PageAllocator(int pgsize_) : pgsize(pgsize_) {
-        new_page();
+    void initialize(int nlocal, int pgsize)
+    {
+        nlocal_ = nlocal;
+        allocator_ = PageAllocator(pgsize);
+        numneigh_.assign(nlocal, 0);
+        firstneigh_.assign(nlocal, nullptr);
     }
 
-    ~PageAllocator() {
-        // no manual delete[]; vectors clean themselves
+    void reset()
+    {
+        allocator_.reset();
     }
 
-    PageAllocator(const PageAllocator&) = delete;
-    PageAllocator& operator=(const PageAllocator&) = delete;
+    int get_nlocal() const { return nlocal_; }
+    int get_numneigh(int i) const { return numneigh_[i]; }
+    int* get_firstneigh(int i) { return firstneigh_[i]; }
+    const int* get_firstneigh(int i) const { return firstneigh_[i]; }
+    PageAllocator& get_allocator() { return allocator_; }
+    const PageAllocator& get_allocator() const { return allocator_; }
 
-    // Allow move
-    PageAllocator(PageAllocator&&) = default;
-    PageAllocator& operator=(PageAllocator&&) = default;
+private:
+    int nlocal_ = 0;
+    std::vector<int> numneigh_;
+    std::vector<int*> firstneigh_;
+    PageAllocator allocator_;
 
-    void new_page() {
-        Page p;
-        p.capacity = pgsize;
-        p.offset = 0;
-        p.data.resize(pgsize);
-        pages.push_back(std::move(p));
-    }
-
-    int* allocate(int n) {
-        if (n <= 0) return nullptr;
-        // reject requests larger than a single page
-        if (n > pgsize) {
-            std::cerr << "PageAllocator::allocate error: request " << n << " larger than page size " << pgsize << std::endl;
-            return nullptr;
-        }
-        if (pages.empty()) new_page();
-        Page& p = pages.back();
-        if (p.offset + n > p.capacity) {
-            new_page();
-            return allocate(n);
-        }
-        int* ptr = p.data.data() + p.offset;
-        p.offset += n;
-        return ptr;
-    }
-
-    void reset() {
-        pages.resize(1);
-        pages[0].offset = 0;
-    }
+    friend class BinManager;
 };
 
-//////////////////////////////////////////////////////////////
-// Neighbor List
-//////////////////////////////////////////////////////////////
-
-class NeighborList {
-public:
-    NeighborList() = default;
-    ~NeighborList() = default;
-
-    int nlocal;
-
-    std::vector<int> numneigh;
-    std::vector<int*> firstneigh;
-
-    PageAllocator allocator;
-
-    void initialize(int n, int pgsize) {
-        nlocal = n;
-    allocator = PageAllocator(pgsize);
-    // ensure neighbor containers are sized and initialized
-    numneigh.assign(n, 0);
-    firstneigh.assign(n, nullptr);
-    }
-
-    void reset() {
-        allocator.reset();
-    }
-};
\ No newline at end of file
+#endif // NEIGHBOR_LIST_H
\ No newline at end of file
diff --git a/source/source_cell/module_neighlist/neighbor_search.cpp b/source/source_cell/module_neighlist/neighbor_search.cpp
index d8f86ecf04c..912515bf9d5 100644
--- a/source/source_cell/module_neighlist/neighbor_search.cpp
+++ b/source/source_cell/module_neighlist/neighbor_search.cpp
@@ -2,9 +2,112 @@
 #include <cmath>
 #include <algorithm>
 #include <limits>
+#include <cassert>
 
+// ========== Getter methods ==========
 
-InputAtoms NeighborSearch::ucell_to_input_atoms(const IAtomProvider& ucell)
+double NeighborSearch::get_search_radius() const {
+    return search_radius_;
+}
+
+int NeighborSearch::get_x() const {
+    return x_;
+}
+
+int NeighborSearch::get_y() const {
+    return y_;
+}
+
+int NeighborSearch::get_z() const {
+    return z_;
+}
+
+double NeighborSearch::get_wide_x() const {
+    return wide_x_;
+}
+
+double NeighborSearch::get_wide_y() const {
+    return wide_y_;
+}
+
+double NeighborSearch::get_wide_z() const {
+    return wide_z_;
+}
+
+int NeighborSearch::get_glayerX() const {
+    return glayerX_;
+}
+
+int NeighborSearch::get_glayerY() const {
+    return glayerY_;
+}
+
+int NeighborSearch::get_glayerZ() const {
+    return glayerZ_;
+}
+
+int NeighborSearch::get_glayerX_minus() const {
+    return glayerX_minus_;
+}
+
+int NeighborSearch::get_glayerY_minus() const {
+    return glayerY_minus_;
+}
+
+int NeighborSearch::get_glayerZ_minus() const {
+    return glayerZ_minus_;
+}
+
+const std::vector<NeighborAtom>& NeighborSearch::get_all_atoms() const {
+    return all_atoms_;
+}
+
+const std::vector<NeighborAtom>& NeighborSearch::get_inside_atoms() const {
+    return inside_atoms_;
+}
+
+const std::vector<NeighborAtom>& NeighborSearch::get_ghost_atoms() const {
+    return ghost_atoms_;
+}
+
+NeighborList& NeighborSearch::get_neighbor_list() {
+    return neighbor_list_;
+}
+
+const NeighborList& NeighborSearch::get_neighbor_list() const {
+    return neighbor_list_;
+}
+
+// ========== Setter methods ==========
+
+void NeighborSearch::set_search_radius(double sr) {
+    search_radius_ = sr;
+}
+
+void NeighborSearch::set_position(int x, int y, int z) {
+    x_ = x;
+    y_ = y;
+    z_ = z;
+}
+
+void NeighborSearch::set_width(double wx, double wy, double wz) {
+    wide_x_ = wx;
+    wide_y_ = wy;
+    wide_z_ = wz;
+}
+
+// ========== Internal methods ==========
+
+double NeighborSearch::cross_product_norm(double a1, double a2, double a3,
+                                          double b1, double b2, double b3)
+{
+    double c1 = a2 * b3 - a3 * b2;
+    double c2 = a3 * b1 - a1 * b3;
+    double c3 = a1 * b2 - a2 * b1;
+    return sqrt(c1 * c1 + c2 * c2 + c3 * c3);
+}
+
+InputAtoms NeighborSearch::ucell_to_input_atoms(const AtomProvider& ucell)
 {
     InputAtoms input_atoms;
     int atom_count = 0;
@@ -18,9 +121,9 @@ InputAtoms NeighborSearch::ucell_to_input_atoms(const IAtomProvider& ucell)
         for (int j = 0; j < ucell.get_na(i); j++)
         {
             NeighborAtom atom(
-                ucell.get_tauu(i,j).x,
-                ucell.get_tauu(i,j).y,
-                ucell.get_tauu(i,j).z,
+                ucell.get_tau(i,j).x,
+                ucell.get_tau(i,j).y,
+                ucell.get_tau(i,j).z,
                 i,
                 j,
                 atom_count
@@ -42,44 +145,112 @@ InputAtoms NeighborSearch::ucell_to_input_atoms(const IAtomProvider& ucell)
     return input_atoms;
 }
 
-void NeighborSearch::init(const IAtomProvider& ucell, double sr, int mpi_rank)
+void NeighborSearch::check_expand_condition(const AtomProvider& ucell)
+{
+    const auto& lat = ucell.get_latvec();
+    const double omega = ucell.get_omega();
+    const double lat0 = ucell.get_lat0();
+    const double lat0_cubed = lat0 * lat0 * lat0;
+
+    double a23_norm = cross_product_norm(lat.e21, lat.e22, lat.e23, lat.e31, lat.e32, lat.e33);
+    int extend_d11 = std::ceil(a23_norm * search_radius_ / omega * lat0_cubed);
+
+    double a31_norm = cross_product_norm(lat.e31, lat.e32, lat.e33, lat.e11, lat.e12, lat.e13);
+    int extend_d22 = std::ceil(a31_norm * search_radius_ / omega * lat0_cubed);
+
+    double a12_norm = cross_product_norm(lat.e11, lat.e12, lat.e13, lat.e21, lat.e22, lat.e23);
+    int extend_d33 = std::ceil(a12_norm * search_radius_ / omega * lat0_cubed);
+
+    glayerX_ = extend_d11 + positive_layer_offset;
+    glayerY_ = extend_d22 + positive_layer_offset;
+    glayerZ_ = extend_d33 + positive_layer_offset;
+    glayerX_minus_ = extend_d11;
+    glayerY_minus_ = extend_d22;
+    glayerZ_minus_ = extend_d33;
+}
+
+void NeighborSearch::set_member_variables(const AtomProvider& ucell)
+{
+    all_atoms_.clear();
+
+    ModuleBase::Vector3<double> vec1(ucell.get_latvec().e11, ucell.get_latvec().e12, ucell.get_latvec().e13);
+    ModuleBase::Vector3<double> vec2(ucell.get_latvec().e21, ucell.get_latvec().e22, ucell.get_latvec().e23);
+    ModuleBase::Vector3<double> vec3(ucell.get_latvec().e31, ucell.get_latvec().e32, ucell.get_latvec().e33);
+
+    int atom_count = 0;
+
+    for (int ix = -glayerX_minus_; ix < glayerX_; ix++)
+    {
+        for (int iy = -glayerY_minus_; iy < glayerY_; iy++)
+        {
+            for (int iz = -glayerZ_minus_; iz < glayerZ_; iz++)
+            {
+                for (int i = 0; i < ucell.get_ntype(); i++)
+                {
+                    for (int j = 0; j < ucell.get_na(i); j++)
+                    {
+                        double atom_x = ucell.get_tau(i,j).x + vec1[0] * ix + vec2[0] * iy + vec3[0] * iz;
+                        double atom_y = ucell.get_tau(i,j).y + vec1[1] * ix + vec2[1] * iy + vec3[1] * iz;
+                        double atom_z = ucell.get_tau(i,j).z + vec1[2] * ix + vec2[2] * iy + vec3[2] * iz;
+
+                        NeighborAtom atom(atom_x, atom_y, atom_z, i, j, atom_count);
+                        if(ix==0 && iy==0 && iz==0)
+                        {
+                            atom.is_inside = true;
+                        }
+                        else
+                        {
+                            atom.is_inside = false;
+                        }
+                        all_atoms_.push_back(atom);
+                        atom_count++;
+                    }
+                }
+            }
+        }
+    }
+}
+
+// ========== Main public interface ==========
+
+void NeighborSearch::init(const AtomProvider& ucell, double sr, int mpi_rank)
 {
     // clear possible residual data from previous runs
-    inside_atoms.clear();
-    ghost_atoms.clear();
-    all_atoms.clear();
+    inside_atoms_.clear();
+    ghost_atoms_.clear();
+    all_atoms_.clear();
     // clear any existing bin manager state
-    bin_manager.clear();
+    bin_manager_.clear();
 
-    search_radius = sr / ucell.get_lat0();
-    Check_Expand_Condition(ucell);
-    setMemberVariables(ucell);
+    search_radius_ = sr / ucell.get_lat0();
+    check_expand_condition(ucell);
+    set_member_variables(ucell);
     InputAtoms atoms = ucell_to_input_atoms(ucell);
 
     int mpi_size = 1;
     int nx, ny, nz;
     decompose(mpi_size, nx, ny, nz);
 
-    z = mpi_rank / (nx * ny);
-    y = (mpi_rank % (nx * ny)) / nx;
-    x = mpi_rank % (nx * ny) % nx;
+    z_ = mpi_rank / (nx * ny);
+    y_ = (mpi_rank % (nx * ny)) / nx;
+    x_ = mpi_rank % (nx * ny) % nx;
 
-    wide_x = (atoms.x_high - atoms.x_low) / nx;
-    wide_y = (atoms.y_high - atoms.y_low) / ny;
-    wide_z = (atoms.z_high - atoms.z_low) / nz;
-    assert(wide_x>=0);
-    assert(wide_y>=0);
-    assert(wide_z>=0);
+    wide_x_ = (atoms.x_high - atoms.x_low) / nx;
+    wide_y_ = (atoms.y_high - atoms.y_low) / ny;
+    wide_z_ = (atoms.z_high - atoms.z_low) / nz;
+    assert(wide_x_ >= 0);
+    assert(wide_y_ >= 0);
+    assert(wide_z_ >= 0);
 
     int in_x, in_y, in_z;
 
-    for (int i = 0; i < all_atoms.size(); i++)
+    for (size_t i = 0; i < all_atoms_.size(); i++)
     {
-        if(wide_x<1e-8)
+        if(wide_x_ < coord_tolerance)
         {
-            if(std::abs(all_atoms[i].position_x-atoms.x_low)<1e-8)
+            if(std::abs(all_atoms_[i].position_x - atoms.x_low) < coord_tolerance)
             {
-                in_x = x;
+                in_x = x_;
             }
             else
             {
@@ -89,15 +260,15 @@ void NeighborSearch::init(const IAtomProvider& ucell, double sr, int mpi_rank)
         else
         {
             in_x = std::min(
-                static_cast<int>(std::floor((all_atoms[i].position_x - atoms.x_low) / wide_x)),
+                static_cast<int>(std::floor((all_atoms_[i].position_x - atoms.x_low) / wide_x_)),
                 nx - 1
             );
         }
-        if(wide_y<1e-8)
+        if(wide_y_ < coord_tolerance)
         {
-            if(std::abs(all_atoms[i].position_y-atoms.y_low)<1e-8)
+            if(std::abs(all_atoms_[i].position_y - atoms.y_low) < coord_tolerance)
             {
-                in_y = y;
+                in_y = y_;
             }
             else
             {
@@ -107,15 +278,15 @@ void NeighborSearch::init(const IAtomProvider& ucell, double sr, int mpi_rank)
         else
         {
             in_y = std::min(
-                static_cast<int>(std::floor((all_atoms[i].position_y - atoms.y_low) / wide_y)),
+                static_cast<int>(std::floor((all_atoms_[i].position_y - atoms.y_low) / wide_y_)),
                 ny - 1
             );
         }
-        if(wide_z<1e-8)
+        if(wide_z_ < coord_tolerance)
         {
-            if(std::abs(all_atoms[i].position_z-atoms.z_low)<1e-8)
+            if(std::abs(all_atoms_[i].position_z - atoms.z_low) < coord_tolerance)
             {
-                in_z = z;
+                in_z = z_;
             }
             else
             {
@@ -125,113 +296,42 @@ void NeighborSearch::init(const IAtomProvider& ucell, double sr, int mpi_rank)
         else
         {
             in_z = std::min(
-                static_cast<int>(std::floor((all_atoms[i].position_z - atoms.z_low) / wide_z)),
+                static_cast<int>(std::floor((all_atoms_[i].position_z - atoms.z_low) / wide_z_)),
                 nz - 1
             );
         }
-        //std::cout<<in_x<<" "<<in_y<<" "<<in_z<<std::endl;
 
-        if (in_x==x && in_y==y && in_z==z&&all_atoms[i].position_x<=atoms.x_high&&all_atoms[i].position_y<=atoms.y_high&&all_atoms[i].position_z<=atoms.z_high&&all_atoms[i].is_inside)
+        if (in_x == x_ && in_y == y_ && in_z == z_ &&
+            all_atoms_[i].position_x <= atoms.x_high &&
+            all_atoms_[i].position_y <= atoms.y_high &&
+            all_atoms_[i].position_z <= atoms.z_high &&
+            all_atoms_[i].is_inside)
         {
-            //all_atoms[i].isghost = false;
-            inside_atoms.push_back(all_atoms[i]);
+            inside_atoms_.push_back(all_atoms_[i]);
         }
         else if (distance(
-            all_atoms[i].position_x,
-            all_atoms[i].position_y,
-            all_atoms[i].position_z,
+            all_atoms_[i].position_x,
+            all_atoms_[i].position_y,
+            all_atoms_[i].position_z,
             atoms.x_low,
             atoms.y_low,
-            atoms.z_low) <= search_radius * search_radius)
+            atoms.z_low) <= search_radius_ * search_radius_)
         {
-            //all_atoms[i].isghost = true;
-            ghost_atoms.push_back(all_atoms[i]);
+            ghost_atoms_.push_back(all_atoms_[i]);
         }
     }
 
-    neighbor_list.initialize(inside_atoms.size(), all_atoms.size()*2);
+    neighbor_list_.initialize(inside_atoms_.size(), all_atoms_.size() * neighbor_reserve_factor);
 }
 
 void NeighborSearch::build_neighbors()
 {
-    bin_manager.init_bins(search_radius, inside_atoms, ghost_atoms);
-    bin_manager.do_binning(inside_atoms, ghost_atoms);
-    bin_manager.build_atom_neighbors(neighbor_list, inside_atoms);
+    bin_manager_.init_bins(search_radius_, inside_atoms_, ghost_atoms_);
+    bin_manager_.do_binning(inside_atoms_, ghost_atoms_);
+    bin_manager_.build_atom_neighbors(neighbor_list_, inside_atoms_);
 }
 
-void NeighborSearch::Check_Expand_Condition(const IAtomProvider& ucell)
-{
-    double a23_1 = ucell.get_latvec().e22 * ucell.get_latvec().e33 - ucell.get_latvec().e23 * ucell.get_latvec().e32;
-    double a23_2 = ucell.get_latvec().e21 * ucell.get_latvec().e33 - ucell.get_latvec().e23 * ucell.get_latvec().e31;
-    double a23_3 = ucell.get_latvec().e21 * ucell.get_latvec().e32 - ucell.get_latvec().e22 * ucell.get_latvec().e31;
-    double a23_norm = sqrt(a23_1 * a23_1 + a23_2 * a23_2 + a23_3 * a23_3);
-    double extend_v = a23_norm * search_radius;
-    double extend_d1 = extend_v / ucell.get_omega() * ucell.get_lat0() * ucell.get_lat0() * ucell.get_lat0();
-    int extend_d11 = std::ceil(extend_d1);
-
-    double a31_1 = ucell.get_latvec().e32 * ucell.get_latvec().e13 - ucell.get_latvec().e33 * ucell.get_latvec().e12;
-    double a31_2 = ucell.get_latvec().e31 * ucell.get_latvec().e13 - ucell.get_latvec().e33 * ucell.get_latvec().e11;
-    double a31_3 = ucell.get_latvec().e31 * ucell.get_latvec().e12 - ucell.get_latvec().e32 * ucell.get_latvec().e11;
-    double a31_norm = sqrt(a31_1 * a31_1 + a31_2 * a31_2 + a31_3 * a31_3);
-    double extend_d2 = a31_norm * search_radius / ucell.get_omega() * ucell.get_lat0() * ucell.get_lat0() * ucell.get_lat0();
-    int extend_d22 = std::ceil(extend_d2);
-
-    double a12_1 = ucell.get_latvec().e12 * ucell.get_latvec().e23 - ucell.get_latvec().e13 * ucell.get_latvec().e22;
-    double a12_2 = ucell.get_latvec().e11 * ucell.get_latvec().e23 - ucell.get_latvec().e13 * ucell.get_latvec().e21;
-    double a12_3 = ucell.get_latvec().e11 * ucell.get_latvec().e22 - ucell.get_latvec().e12 * ucell.get_latvec().e21;
-    double a12_norm = sqrt(a12_1 * a12_1 + a12_2 * a12_2 + a12_3 * a12_3);
-    double extend_d3 = a12_norm * search_radius / ucell.get_omega() * ucell.get_lat0() * ucell.get_lat0() * ucell.get_lat0();
-    int extend_d33 = std::ceil(extend_d3);
-
-    glayerX = extend_d11 + 1;
-    glayerY = extend_d22 + 1;
-    glayerZ = extend_d33 + 1;
-    glayerX_minus = extend_d11;
-    glayerY_minus = extend_d22;
-    glayerZ_minus = extend_d33;
-}
-
-void NeighborSearch::setMemberVariables(const IAtomProvider& ucell)
-{
-    all_atoms.clear();
-
-    ModuleBase::Vector3<double> vec1(ucell.get_latvec().e11, ucell.get_latvec().e12, ucell.get_latvec().e13);
-    ModuleBase::Vector3<double> vec2(ucell.get_latvec().e21, ucell.get_latvec().e22, ucell.get_latvec().e23);
-    ModuleBase::Vector3<double> vec3(ucell.get_latvec().e31, ucell.get_latvec().e32, ucell.get_latvec().e33);
-
-    int atom_count = 0;
-
-    for (int ix = -glayerX_minus; ix < glayerX; ix++)
-    {
-        for (int iy = -glayerY_minus; iy < glayerY; iy++)
-        {
-            for (int iz = -glayerZ_minus; iz < glayerZ; iz++)
-            {
-                for (int i = 0; i < ucell.get_ntype(); i++)
-                {
-                    for (int j = 0; j < ucell.get_na(i); j++)
-                    {
-                        double x = ucell.get_tauu(i,j).x + vec1[0] * ix + vec2[0] * iy + vec3[0] * iz;
-                        double y = ucell.get_tauu(i,j).y + vec1[1] * ix + vec2[1] * iy + vec3[1] * iz;
-                        double z = ucell.get_tauu(i,j).z + vec1[2] * ix + vec2[2] * iy + vec3[2] * iz;
-
-                        NeighborAtom atom(x, y, z, i, j, atom_count);
-                        if(ix==0&&iy==0&&iz==0)
-                        {
-                            atom.is_inside = true;
-                        }
-                        else
-                        {
-                            atom.is_inside = false;
-                        }
-                        all_atoms.push_back(atom);
-                        atom_count++;
-                    }
-                }
-            }
-        }
-    }
-}
+// ========== Utility methods ==========
 
 double NeighborSearch::distance(
     double position_x,
@@ -241,9 +341,9 @@ double NeighborSearch::distance(
     double y_low,
     double z_low)
 {
-    double dx = std::max(0.0, std::max(x_low + x * wide_x - position_x, position_x - (x_low + (x + 1) * wide_x)));
-    double dy = std::max(0.0, std::max(y_low + y * wide_y - position_y, position_y - (y_low + (y + 1) * wide_y)));
-    double dz = std::max(0.0, std::max(z_low + z * wide_z - position_z, position_z - (z_low + (z + 1) * wide_z)));
+    double dx = std::max(0.0, std::max(x_low + x_ * wide_x_ - position_x, position_x - (x_low + (x_ + 1) * wide_x_)));
+    double dy = std::max(0.0, std::max(y_low + y_ * wide_y_ - position_y, position_y - (y_low + (y_ + 1) * wide_y_)));
+    double dz = std::max(0.0, std::max(z_low + z_ * wide_z_ - position_z, position_z - (z_low + (z_ + 1) * wide_z_)));
     return dx * dx + dy * dy + dz * dz;
 }
 
@@ -253,7 +353,7 @@ void NeighborSearch::decompose(int mpi_size, int &nx, int &ny, int &nz)
     ny = 1;
     nz = mpi_size;
 
-    int cube = cbrt(mpi_size);
+    int cube = static_cast<int>(cbrt(mpi_size));
     for (int i = cube; i >= 1; i--)
     {
         if (mpi_size % i == 0)
@@ -264,7 +364,7 @@ void NeighborSearch::decompose(int mpi_size, int &nx, int &ny, int &nz)
         }
     }
 
-    int sq = sqrt(ny);
+    int sq = static_cast<int>(sqrt(ny));
     for (int i = sq; i >= 1; i--)
     {
         if (ny % i == 0)
diff --git a/source/source_cell/module_neighlist/neighbor_search.h b/source/source_cell/module_neighlist/neighbor_search.h
index 53b36761015..b75a8926d5c 100644
--- a/source/source_cell/module_neighlist/neighbor_search.h
+++ b/source/source_cell/module_neighlist/neighbor_search.h
@@ -1,55 +1,316 @@
 #ifndef NEIGHBOR_SEARCH_H
 #define NEIGHBOR_SEARCH_H
 
-
 #include "source_cell/module_neighlist/neighbor_atom.h"
 #include "source_cell/module_neighlist/bin_manager.h"
 #include "source_cell/module_neighlist/neighbor_list.h"
-#include "source_cell/module_neighlist/unitcell_interface.h"
+#include "source_cell/module_neighlist/atom_provider.h"
 
+/**
+ * @brief Neighbor search algorithm for building atom neighbor lists.
+ *
+ * This class implements a neighbor search algorithm that finds all atoms
+ * within a given cutoff radius. It uses a binning strategy for efficiency
+ * and supports MPI parallelization by decomposing the simulation domain.
+ *
+ * The workflow is:
+ * 1. Call init() to initialize with unit cell and search radius
+ * 2. Call build_neighbors() to construct the neighbor list
+ * 3. Access results via get_neighbor_list()
+ */
 class NeighborSearch
 {
 public:
-    NeighborSearch()=default;
-    ~NeighborSearch()=default;
+    /**
+     * @brief Default constructor.
+     */
+    NeighborSearch() = default;
+
+    /**
+     * @brief Default destructor.
+     */
+    ~NeighborSearch() = default;
 
-    void init(const IAtomProvider& ucell, double sr, int mpi_rank);
+    // ========== Main public interface ==========
 
+    /**
+     * @brief Initialize the neighbor search with unit cell and search radius.
+     *
+     * This method sets up the domain decomposition, identifies inside and
+     * ghost atoms, and prepares internal data structures.
+     *
+     * @param ucell Unit cell providing atom positions and lattice info.
+     * @param sr Search radius (cutoff distance) in Bohr.
+     * @param mpi_rank MPI rank of this process.
+     */
+    void init(const AtomProvider& ucell, double sr, int mpi_rank);
 
+    /**
+     * @brief Build the neighbor list for all inside atoms.
+     *
+     * Must be called after init(). Uses binning to efficiently find
+     * all neighbors within the search radius.
+     */
     void build_neighbors();
 
-    InputAtoms ucell_to_input_atoms(const IAtomProvider& ucell);
+    /**
+     * @brief Get the constructed neighbor list.
+     * @return Reference to the NeighborList object.
+     */
+    NeighborList& get_neighbor_list();
 
-    void Check_Expand_Condition(const IAtomProvider& ucell);
+    /**
+     * @brief Get the constructed neighbor list (const version).
+     * @return Const reference to the NeighborList object.
+     */
+    const NeighborList& get_neighbor_list() const;
 
-    void setMemberVariables(const IAtomProvider& ucell);
-    NeighborList& get_neighbor_list() { return neighbor_list; }
+    // ========== Utility methods (public for testing) ==========
 
-    double distance(
-        double position_x,
-        double position_y,
-        double position_z,
-        double x_low,
-        double y_low,
-        double z_low);
+    /**
+     * @brief Calculate squared distance from a point to the local domain box.
+     *
+     * Used to determine if an atom is within the search radius of the
+     * local MPI domain.
+     *
+     * @param position_x X coordinate of the point.
+     * @param position_y Y coordinate of the point.
+     * @param position_z Z coordinate of the point.
+     * @param x_low Lower bound of the global domain in X.
+     * @param y_low Lower bound of the global domain in Y.
+     * @param z_low Lower bound of the global domain in Z.
+     * @return Squared distance to the domain box.
+     */
+    double distance(double position_x,
+                    double position_y,
+                    double position_z,
+                    double x_low,
+                    double y_low,
+                    double z_low);
 
+    /**
+     * @brief Decompose MPI size into a 3D grid.
+     *
+     * Finds a balanced decomposition of mpi_size into nx * ny * nz.
+     *
+     * @param mpi_size Total number of MPI processes.
+     * @param nx Output: number of divisions in X.
+     * @param ny Output: number of divisions in Y.
+     * @param nz Output: number of divisions in Z.
+     */
     void decompose(int mpi_size, int& nx, int& ny, int& nz);
 
-    double search_radius;
+    // ========== Getter methods ==========
+
+    /**
+     * @brief Get the search radius.
+     * @return Search radius in lattice units.
+     */
+    double get_search_radius() const;
+
+    /**
+     * @brief Get the X position of this MPI domain.
+     * @return Domain index in X.
+     */
+    int get_x() const;
+
+    /**
+     * @brief Get the Y position of this MPI domain.
+     * @return Domain index in Y.
+     */
+    int get_y() const;
+
+    /**
+     * @brief Get the Z position of this MPI domain.
+     * @return Domain index in Z.
+     */
+    int get_z() const;
+
+    /**
+     * @brief Get the width of this MPI domain in X.
+     * @return Domain width in X.
+     */
+    double get_wide_x() const;
+
+    /**
+     * @brief Get the width of this MPI domain in Y.
+     * @return Domain width in Y.
+     */
+    double get_wide_y() const;
+
+    /**
+     * @brief Get the width of this MPI domain in Z.
+     * @return Domain width in Z.
+     */
+    double get_wide_z() const;
+
+    /**
+     * @brief Get the number of expansion layers in +X direction.
+     * @return Number of layers.
+     */
+    int get_glayerX() const;
+
+    /**
+     * @brief Get the number of expansion layers in +Y direction.
+     * @return Number of layers.
+     */
+    int get_glayerY() const;
+
+    /**
+     * @brief Get the number of expansion layers in +Z direction.
+     * @return Number of layers.
+     */
+    int get_glayerZ() const;
+
+    /**
+     * @brief Get the number of expansion layers in -X direction.
+     * @return Number of layers.
+     */
+    int get_glayerX_minus() const;
+
+    /**
+     * @brief Get the number of expansion layers in -Y direction.
+     * @return Number of layers.
+     */
+    int get_glayerY_minus() const;
+
+    /**
+     * @brief Get the number of expansion layers in -Z direction.
+     * @return Number of layers.
+     */
+    int get_glayerZ_minus() const;
+
+    /**
+     * @brief Get all atoms (including periodic images).
+     * @return Const reference to the vector of all atoms.
+     */
+    const std::vector<NeighborAtom>& get_all_atoms() const;
+
+    /**
+     * @brief Get atoms inside the local MPI domain.
+     * @return Const reference to the vector of inside atoms.
+     */
+    const std::vector<NeighborAtom>& get_inside_atoms() const;
+
+    /**
+     * @brief Get ghost atoms (neighbors of inside atoms).
+     * @return Const reference to the vector of ghost atoms.
+     */
+    const std::vector<NeighborAtom>& get_ghost_atoms() const;
+
+    // ========== Setter methods ==========
+
+    /**
+     * @brief Set the search radius.
+     * @param sr Search radius in lattice units.
+     */
+    void set_search_radius(double sr);
+
+    /**
+     * @brief Set the position of this MPI domain.
+     * @param x Domain index in X.
+     * @param y Domain index in Y.
+     * @param z Domain index in Z.
+     */
+    void set_position(int x, int y, int z);
+
+    /**
+     * @brief Set the width of this MPI domain.
+     * @param wx Domain width in X.
+     * @param wy Domain width in Y.
+     * @param wz Domain width in Z.
+     */
+    void set_width(double wx, double wy, double wz);
+
+private:
+    // ========== Internal methods ==========
+
+    /**
+     * @brief Convert unit cell atoms to InputAtoms format.
+     * @param ucell Unit cell providing atom info.
+     * @return InputAtoms structure for processing.
+     */
+    InputAtoms ucell_to_input_atoms(const AtomProvider& ucell);
+
+    /**
+     * @brief Check and compute expansion layer counts.
+     *
+     * Determines how many periodic images are needed to cover
+     * the search radius in each lattice direction.
+     *
+     * @param ucell Unit cell providing lattice vectors.
+     */
+    void check_expand_condition(const AtomProvider& ucell);
+
+    /**
+     * @brief Set member variables by generating periodic images.
+     *
+     * Populates all_atoms_ with atoms from the unit cell and
+     * all required periodic images.
+     *
+     * @param ucell Unit cell providing atom positions.
+     */
+    void set_member_variables(const AtomProvider& ucell);
+
+    /**
+     * @brief Compute the norm of the cross product of two 3D vectors.
+     *
+     * @param a1, a2, a3 Components of the first vector.
+     * @param b1, b2, b3 Components of the second vector.
+     * @return Norm of the cross product.
+     */
+    static double cross_product_norm(double a1, double a2, double a3,
+                                     double b1, double b2, double b3);
+
+    // ========== Data members ==========
+
+    /// Search radius in lattice units
+    double search_radius_ = 0.0;
+
+    /// Position of this MPI domain in the 3D grid
+    int x_ = 0;
+    int y_ = 0;
+    int z_ = 0;
+
+    /// Width of this MPI domain
+    double wide_x_ = 0.0;
+    double wide_y_ = 0.0;
+    double wide_z_ = 0.0;
+
+    /// Number of expansion layers in positive directions
+    int glayerX_ = 0;
+    int glayerY_ = 0;
+    int glayerZ_ = 0;
+
+    /// Number of expansion layers in negative directions
+    int glayerX_minus_ = 0;
+    int glayerY_minus_ = 0;
+    int glayerZ_minus_ = 0;
+
+    /// All atoms including periodic images
+    std::vector<NeighborAtom> all_atoms_;
+
+    /// Atoms inside the local MPI domain
+    std::vector<NeighborAtom> inside_atoms_;
+
+    /// Ghost atoms (neighbors from other domains or images)
+    std::vector<NeighborAtom> ghost_atoms_;
 
+    /// The constructed neighbor list
+    NeighborList neighbor_list_;
 
-    int x, y, z;
-    double wide_x, wide_y, wide_z;
+    /// Bin manager for efficient neighbor search
+    BinManager bin_manager_;
+    // ========== Compile-time constants ==========
 
-    int glayerX, glayerY, glayerZ;
-    int glayerX_minus, glayerY_minus, glayerZ_minus;
+    /// Tolerance for coordinate comparisons in lattice units
+    static constexpr double coord_tolerance = 1e-8;
 
-    std::vector<NeighborAtom> all_atoms;
-    std::vector<NeighborAtom> inside_atoms;
-    std::vector<NeighborAtom> ghost_atoms;
+    /// Offset added to expansion layers in positive directions
+    static constexpr int positive_layer_offset = 1;
 
-    NeighborList neighbor_list;
-    BinManager bin_manager;
+    /// Reserve factor for neighbor list capacity estimation
+    static constexpr int neighbor_reserve_factor = 2;
 };
 
-#endif
\ No newline at end of file
+#endif // NEIGHBOR_SEARCH_H
\ No newline at end of file
diff --git a/source/source_cell/module_neighlist/page_allocator.cpp b/source/source_cell/module_neighlist/page_allocator.cpp
new file mode 100644
index 00000000000..959ea79154f
--- /dev/null
+++ b/source/source_cell/module_neighlist/page_allocator.cpp
@@ -0,0 +1,68 @@
+#include "page_allocator.h"
+#include "source_base/tool_quit.h"
+
+PageAllocator::PageAllocator() : pgsize_(default_pgsize)
+{
+    new_page_();
+}
+
+PageAllocator::PageAllocator(int pgsize) : pgsize_(pgsize)
+{
+    new_page_();
+}
+
+PageAllocator::~PageAllocator() = default;
+
+int* PageAllocator::allocate(int n)
+{
+    if (n <= 0) 
+    {
+        return nullptr;
+    }
+
+    if (n > pgsize_)
+    {
+        ModuleBase::WARNING_QUIT(
+            "PageAllocator::allocate",
+            "request " + std::to_string(n) + " larger than page size " + std::to_string(pgsize_)
+        );
+    }
+
+    if (pages_.empty()) 
+    {
+        new_page_();
+    }
+
+    Page& p = pages_.back();
+
+    if (p.offset + n > p.capacity)
+    {
+        new_page_();
+        return allocate(n);
+    }
+
+    int* ptr = p.data.data() + p.offset;
+    p.offset += n;
+
+    return ptr;
+}
+
+void PageAllocator::reset()
+{
+    pages_.resize(1);
+    pages_[0].offset = 0;
+}
+
+int PageAllocator::get_pgsize() const
+{
+    return pgsize_;
+}
+
+void PageAllocator::new_page_()
+{
+    Page p;
+    p.capacity = pgsize_;
+    p.offset = 0;
+    p.data.resize(pgsize_);
+    pages_.push_back(std::move(p));
+}
diff --git a/source/source_cell/module_neighlist/page_allocator.h b/source/source_cell/module_neighlist/page_allocator.h
new file mode 100644
index 00000000000..e6cb9e6756f
--- /dev/null
+++ b/source/source_cell/module_neighlist/page_allocator.h
@@ -0,0 +1,38 @@
+#ifndef PAGE_ALLOCATOR_H
+#define PAGE_ALLOCATOR_H
+
+#include <vector>
+
+class PageAllocator
+{
+public:
+    enum { default_pgsize = 1024 };
+
+    PageAllocator();
+    explicit PageAllocator(int pgsize);
+    ~PageAllocator();
+
+    PageAllocator(const PageAllocator&) = delete;
+    PageAllocator& operator=(const PageAllocator&) = delete;
+    PageAllocator(PageAllocator&&) = default;
+    PageAllocator& operator=(PageAllocator&&) = default;
+
+    int* allocate(int n);
+    void reset();
+    int get_pgsize() const;
+
+private:
+    struct Page
+    {
+        int capacity = 0;
+        int offset = 0;
+        std::vector<int> data;
+    };
+
+    std::vector<Page> pages_;
+    int pgsize_ = 0;
+
+    void new_page_();
+};
+
+#endif // PAGE_ALLOCATOR_H
\ No newline at end of file
diff --git a/source/source_cell/module_neighlist/test/CMakeLists.txt b/source/source_cell/module_neighlist/test/CMakeLists.txt
index a8462c59f7f..ec6df835d5e 100644
--- a/source/source_cell/module_neighlist/test/CMakeLists.txt
+++ b/source/source_cell/module_neighlist/test/CMakeLists.txt
@@ -13,6 +13,8 @@ AddTest(
     neighbor_search_test.cpp
   ../neighbor_search.cpp
   ../bin_manager.cpp
+  ../page_allocator.cpp
+  ../unitcell_lite.cpp
 )
 
 AddTest(
@@ -21,6 +23,7 @@ AddTest(
   SOURCES
     bin_manager_test.cpp
     ../bin_manager.cpp
+    ../page_allocator.cpp
 )
 
 AddTest(
@@ -28,6 +31,7 @@ AddTest(
   LIBS parameter ${math_libs} base device
   SOURCES
   neighbor_list_test.cpp
+  ../page_allocator.cpp
 )
 
 
diff --git a/source/source_cell/module_neighlist/test/bin_manager_test.cpp b/source/source_cell/module_neighlist/test/bin_manager_test.cpp
index 51e5e5e7400..d28b386afb1 100644
--- a/source/source_cell/module_neighlist/test/bin_manager_test.cpp
+++ b/source/source_cell/module_neighlist/test/bin_manager_test.cpp
@@ -7,27 +7,24 @@ TEST(BinManagerUnit, InitAndBinning)
     std::vector<NeighborAtom> inside;
     std::vector<NeighborAtom> ghost;
 
-    // two atoms close to each other
     inside.emplace_back(0.0, 0.0, 0.0, 0, 0, 0);
     inside.emplace_back(0.5, 0.0, 0.0, 0, 1, 1);
 
     BinManager bm;
     bm.init_bins(1.0, inside, ghost);
 
-    EXPECT_EQ(bm.nbinx, 1);
-    EXPECT_EQ(bm.nbiny, 1);
-    EXPECT_EQ(bm.nbinz, 1);
-    EXPECT_EQ(static_cast<int>(bm.bins.size()), bm.nbinx * bm.nbiny * bm.nbinz);
+    EXPECT_EQ(bm.get_nbinx(), 1);
+    EXPECT_EQ(bm.get_nbiny(), 1);
+    EXPECT_EQ(bm.get_nbinz(), 1);
+    EXPECT_EQ(bm.get_total_bins(), bm.get_nbinx() * bm.get_nbiny() * bm.get_nbinz());
 
     bm.do_binning(inside, ghost);
 
-    // compute bin index for first atom
-    int ix = std::min(std::max(int((inside[0].position_x - bm.x_min) / bm.bin_sizex), 0), bm.nbinx - 1);
-    int iy = std::min(std::max(int((inside[0].position_y - bm.y_min) / bm.bin_sizey), 0), bm.nbiny - 1);
-    int iz = std::min(std::max(int((inside[0].position_z - bm.z_min) / bm.bin_sizez), 0), bm.nbinz - 1);
-    int idx = ix * bm.nbiny * bm.nbinz + iy * bm.nbinz + iz;
-
-    EXPECT_GE(bm.bins[idx].atoms.size(), 1u);
+    int total_atoms_in_bins = 0;
+    for (int i = 0; i < bm.get_total_bins(); ++i) {
+        total_atoms_in_bins += bm.get_bin_atom_count(i);
+    }
+    EXPECT_GE(total_atoms_in_bins, 2);
 }
 
 TEST(BinManagerUnit, InitBins)
@@ -42,9 +39,9 @@ TEST(BinManagerUnit, InitBins)
 
     BinManager bm;
     bm.init_bins(1.0, inside, ghost);
-    EXPECT_EQ(bm.nbinx, 5);
-    EXPECT_EQ(bm.nbiny, 1);
-    EXPECT_EQ(bm.nbinz, 1);
+    EXPECT_EQ(bm.get_nbinx(), 5);
+    EXPECT_EQ(bm.get_nbiny(), 1);
+    EXPECT_EQ(bm.get_nbinz(), 1);
 }
 
 TEST(BinManagerUnit, BuildNeighborsAndClear)
@@ -59,10 +56,10 @@ TEST(BinManagerUnit, BuildNeighborsAndClear)
 
     BinManager bm;
     bm.init_bins(1.0, inside, ghost);
-    EXPECT_EQ(bm.nbinx, 5);
-    EXPECT_EQ(bm.nbiny, 1);
-    EXPECT_EQ(bm.nbinz, 1);
-    EXPECT_EQ(static_cast<int>(bm.bins.size()), bm.nbinx * bm.nbiny * bm.nbinz);
+    EXPECT_EQ(bm.get_nbinx(), 5);
+    EXPECT_EQ(bm.get_nbiny(), 1);
+    EXPECT_EQ(bm.get_nbinz(), 1);
+    EXPECT_EQ(bm.get_total_bins(), bm.get_nbinx() * bm.get_nbiny() * bm.get_nbinz());
 
     bm.do_binning(inside, ghost);
 
@@ -71,13 +68,12 @@ TEST(BinManagerUnit, BuildNeighborsAndClear)
 
     bm.build_atom_neighbors(nl, atoms);
 
-    // atom 0 and 1 are neighbors; atom 2 is far
-    EXPECT_EQ(nl.numneigh[0], 1);
-    EXPECT_EQ(nl.numneigh[1], 1);
-    EXPECT_EQ(nl.numneigh[2], 0);
+    EXPECT_EQ(nl.get_numneigh(0), 1);
+    EXPECT_EQ(nl.get_numneigh(1), 1);
+    EXPECT_EQ(nl.get_numneigh(2), 0);
 
     bm.clear();
-    EXPECT_EQ(bm.bins.size(), 0u);
+    EXPECT_EQ(bm.get_total_bins(), 0);
 }
 
 TEST(BinManagerUnit, EmptyAtomsBuildNeighbors)
@@ -91,19 +87,15 @@ TEST(BinManagerUnit, EmptyAtomsBuildNeighbors)
     NeighborList nl;
     nl.initialize(0, 16);
 
-    // should not crash or assert when atoms is empty
     bm.build_atom_neighbors(nl, atoms);
-    EXPECT_EQ(nl.numneigh.size(), 0);
+    EXPECT_EQ(nl.get_nlocal(), 0);
 }
 
 TEST(BinManagerUnit, BoundaryAndExactRadius)
 {
-    // inside atom at origin; other atoms placed on bin boundaries and at exact radius
     std::vector<NeighborAtom> atoms;
     atoms.emplace_back(0.0, 0.0, 0.0, 0, 0, 0);
-    // exactly at search radius 1.0 along x direction
     atoms.emplace_back(1.0, 0.0, 0.0, 0, 1, 1);
-    // slightly inside
     atoms.emplace_back(0.9, 0.0, 0.0, 0, 2, 2);
 
     std::vector<NeighborAtom> inside = atoms;
@@ -118,12 +110,10 @@ TEST(BinManagerUnit, BoundaryAndExactRadius)
 
     bm.build_atom_neighbors(nl, inside);
 
-    // atom at distance exactly 1.0 should be considered neighbor (dist2 == sradius2)
-    EXPECT_EQ(nl.numneigh[0], 2);
-    // the exact atom itself must not be counted as its own neighbor
+    EXPECT_EQ(nl.get_numneigh(0), 2);
     for (int i = 0; i < static_cast<int>(inside.size()); ++i) {
-        for (int j = 0; j < nl.numneigh[i]; ++j) {
-            int id = nl.firstneigh[i][j];
+        for (int j = 0; j < nl.get_numneigh(i); ++j) {
+            int id = nl.get_firstneigh(i)[j];
             EXPECT_NE(id, inside[i].atom_id);
         }
     }
@@ -131,7 +121,6 @@ TEST(BinManagerUnit, BoundaryAndExactRadius)
 
 TEST(BinManagerUnit, InitWithGhostOnly)
 {
-    // inside empty, ghosts present -> init_bins should still compute bounds from ghosts
     std::vector<NeighborAtom> inside;
     std::vector<NeighborAtom> ghost;
 
@@ -141,14 +130,13 @@ TEST(BinManagerUnit, InitWithGhostOnly)
     BinManager bm;
     bm.init_bins(1.0, inside, ghost);
 
-    EXPECT_EQ(bm.nbinx, 3);
-    EXPECT_EQ(bm.nbiny, 1);
-    EXPECT_EQ(bm.nbinz, 1);
+    EXPECT_EQ(bm.get_nbinx(), 3);
+    EXPECT_EQ(bm.get_nbiny(), 1);
+    EXPECT_EQ(bm.get_nbinz(), 1);
 }
 
 TEST(BinManagerUnit, BuildNeighborsNoNeighborsFirstneighNull)
 {
-    // atoms all far apart => no neighbors; firstneigh entries should be nullptr
     std::vector<NeighborAtom> atoms;
     atoms.emplace_back(0.0, 0.0, 0.0, 0, 0, 0);
     atoms.emplace_back(100.0, 100.0, 100.0, 0, 1, 1);
@@ -165,15 +153,14 @@ TEST(BinManagerUnit, BuildNeighborsNoNeighborsFirstneighNull)
 
     bm.build_atom_neighbors(nl, inside);
 
-    EXPECT_EQ(nl.numneigh[0], 0);
-    EXPECT_EQ(nl.numneigh[1], 0);
-    EXPECT_EQ(nl.firstneigh[0], nullptr);
-    EXPECT_EQ(nl.firstneigh[1], nullptr);
+    EXPECT_EQ(nl.get_numneigh(0), 0);
+    EXPECT_EQ(nl.get_numneigh(1), 0);
+    EXPECT_EQ(nl.get_firstneigh(0), nullptr);
+    EXPECT_EQ(nl.get_firstneigh(1), nullptr);
 }
 
 TEST(BinManagerUnit, GhostAtomsAreCounted)
 {
-    // inside atom at origin; ghost atom within search radius should be found
     std::vector<NeighborAtom> inside;
     std::vector<NeighborAtom> ghost;
 
@@ -189,13 +176,12 @@ TEST(BinManagerUnit, GhostAtomsAreCounted)
 
     bm.build_atom_neighbors(nl, inside);
 
-    EXPECT_EQ(nl.numneigh.size(), 1);
-    EXPECT_EQ(nl.numneigh[0], 1);
-    // ensure neighbor id matches ghost atom id
+    EXPECT_EQ(nl.get_nlocal(), 1);
+    EXPECT_EQ(nl.get_numneigh(0), 1);
     bool found = false;
-    if (nl.numneigh[0] > 0 && nl.firstneigh[0] != nullptr) {
-        for (int k = 0; k < nl.numneigh[0]; ++k) {
-            if (nl.firstneigh[0][k] == 3) found = true;
+    if (nl.get_numneigh(0) > 0 && nl.get_firstneigh(0) != nullptr) {
+        for (int k = 0; k < nl.get_numneigh(0); ++k) {
+            if (nl.get_firstneigh(0)[k] == 3) found = true;
         }
     }
     EXPECT_TRUE(found);
@@ -203,7 +189,6 @@ TEST(BinManagerUnit, GhostAtomsAreCounted)
 
 TEST(BinManagerUnit, MultipleBinsNeighborSearch)
 {
-    // Create a 3x3x3 grid of atoms spaced by 1.0 so multiple bins exist
     std::vector<NeighborAtom> atoms;
     int id = 0;
     for (int x = 0; x < 3; ++x)
@@ -223,7 +208,6 @@ TEST(BinManagerUnit, MultipleBinsNeighborSearch)
 
     bm.build_atom_neighbors(nl, inside);
 
-    // center atom (1,1,1) should have multiple neighbors
-    int center_index = 13; // 1*9 + 1*3 + 1
-    EXPECT_EQ(nl.numneigh[center_index], 6);
-}
+    int center_index = 13;
+    EXPECT_EQ(nl.get_numneigh(center_index), 6);
+}
\ No newline at end of file
diff --git a/source/source_cell/module_neighlist/test/neighbor_list_test.cpp b/source/source_cell/module_neighlist/test/neighbor_list_test.cpp
index 8a9b7160543..1731e82c352 100644
--- a/source/source_cell/module_neighlist/test/neighbor_list_test.cpp
+++ b/source/source_cell/module_neighlist/test/neighbor_list_test.cpp
@@ -1,28 +1,20 @@
 #include <gtest/gtest.h>
 #include "../neighbor_list.h"
 
-// Full coverage tests for PageAllocator and NeighborList
-
-TEST(PageAllocator_ConstructorsAndNewPage, DefaultAndCustom)
+TEST(PageAllocator_Constructors, DefaultAndCustom)
 {
-    PageAllocator pa_def; // default page size
-    ASSERT_EQ(pa_def.pages.size(), 1u);
-    EXPECT_EQ(pa_def.pages[0].capacity, pa_def.pgsize);
+    PageAllocator pa_def;
+    EXPECT_EQ(pa_def.get_pgsize(), PageAllocator::default_pgsize);
 
     PageAllocator pa(4);
-    EXPECT_EQ(pa.pgsize, 4);
-    size_t before = pa.pages.size();
-    pa.new_page();
-    EXPECT_EQ(pa.pages.size(), before + 1);
+    EXPECT_EQ(pa.get_pgsize(), 4);
 }
 
-TEST(PageAllocator_AllocateEdgeCases, ZeroNegativeAndTooLarge)
+TEST(PageAllocator_AllocateEdgeCases, ZeroNegative)
 {
     PageAllocator pa(8);
     EXPECT_EQ(pa.allocate(0), nullptr);
     EXPECT_EQ(pa.allocate(-5), nullptr);
-    // larger than single page -> rejected
-    EXPECT_EQ(pa.allocate(9), nullptr);
 }
 
 TEST(PageAllocator_AllocationBehavior, ExactPageAndOverflow)
@@ -42,50 +34,39 @@ TEST(PageAllocator_AllocationBehavior, ExactPageAndOverflow)
     EXPECT_NE(b, a + 2);
 }
 
-TEST(PageAllocator_ResetAndPages, ClearAndReset)
+TEST(PageAllocator_Reset, ClearAndReset)
 {
     PageAllocator pa(4);
     pa.allocate(3);
     pa.allocate(3);
-    ASSERT_EQ(pa.pages.size(), 2u);
 
     pa.reset();
-    EXPECT_EQ(pa.pages.size(), 1u);
-    EXPECT_EQ(pa.pages[0].offset, 0);
-
-    pa.pages.clear();
     int* p = pa.allocate(1);
     ASSERT_NE(p, nullptr);
-    EXPECT_EQ(pa.pages.back().offset, 1);
 }
 
 TEST(NeighborList_InitializeAndReset, Behavior)
 {
     NeighborList nl;
     nl.initialize(0, 16);
-    EXPECT_EQ(nl.nlocal, 0);
-    EXPECT_EQ(nl.numneigh.size(), 0u);
-    EXPECT_EQ(nl.firstneigh.size(), 0u);
+    EXPECT_EQ(nl.get_nlocal(), 0);
 
     nl.initialize(5, 8);
-    EXPECT_EQ(nl.nlocal, 5);
-    EXPECT_EQ(nl.numneigh.size(), 5u);
-    EXPECT_EQ(nl.firstneigh.size(), 5u);
+    EXPECT_EQ(nl.get_nlocal(), 5);
     for (int i = 0; i < 5; ++i) {
-        EXPECT_EQ(nl.numneigh[i], 0);
-        EXPECT_EQ(nl.firstneigh[i], nullptr);
+        EXPECT_EQ(nl.get_numneigh(i), 0);
+        EXPECT_EQ(nl.get_firstneigh(i), nullptr);
     }
-    EXPECT_EQ(nl.allocator.pgsize, 8);
-
-    int* storage = nl.allocator.allocate(3);
-    ASSERT_NE(storage, nullptr);
-    nl.numneigh[2] = 3;
-    nl.firstneigh[2] = storage;
 
     nl.reset();
-    EXPECT_EQ(nl.allocator.pages.size(), 1u);
-    EXPECT_EQ(nl.allocator.pages[0].offset, 0);
-    EXPECT_EQ(nl.numneigh.size(), 5u);
+    EXPECT_EQ(nl.get_nlocal(), 5);
 }
 
-
+TEST(NeighborList_Getters, Accessors)
+{
+    NeighborList nl;
+    nl.initialize(3, 16);
+    EXPECT_EQ(nl.get_nlocal(), 3);
+    EXPECT_EQ(nl.get_numneigh(0), 0);
+    EXPECT_EQ(nl.get_firstneigh(0), nullptr);
+}
\ No newline at end of file
diff --git a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp
index 3387f7fe60c..b8a5bdf0ef3 100644
--- a/source/source_cell/module_neighlist/test/neighbor_search_test.cpp
+++ b/source/source_cell/module_neighlist/test/neighbor_search_test.cpp
@@ -1,30 +1,31 @@
 #include <gtest/gtest.h>
 #include "../neighbor_search.h"
-#include "../unitcell_plus.h"
+#include "../unitcell_lite.h"
+
+// Helper function to create a simple UnitCellLite for testing
+static UnitCellLite make_test_ucell(double lat0, double omega,
+                                     const ModuleBase::Matrix3& latvec,
+                                     int ntype, const std::vector<int>& na,
+                                     const std::vector<ModuleBase::Vector3<double>>& tau) {
+    UnitCellLite ucell;
+    ucell.set_lattice(lat0, omega, latvec);
+    ucell.set_atoms(ntype, na, tau);
+    return ucell;
+}
 
 TEST(NeighborSearchTest, TwoAtomsNeighbor)
 {
-    UnitCellPlus ucell;
-
-    ucell.lat0 = 1.0;
-    ucell.omega = 1.0;
-
-    ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0;
-    ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0;
-    ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1;
+    ModuleBase::Matrix3 latvec;
+    latvec.e11 = 1; latvec.e12 = 0; latvec.e13 = 0;
+    latvec.e21 = 0; latvec.e22 = 1; latvec.e23 = 0;
+    latvec.e31 = 0; latvec.e32 = 0; latvec.e33 = 1;
 
-    ucell.ntype = 1;
-    ucell.na = {2};
-    ucell.nat = 2;
-
-    ucell.tau = {
-        {0.0, 0.0, 0.0},
-        {0.5, 0.0, 0.0}
-    };
-    ucell.compute_naa();
+    UnitCellLite ucell = make_test_ucell(
+        1.0, 1.0, latvec, 1, {2},
+        {{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}
+    );
 
     NeighborSearch ns;
-
     double cutoff = 1.0;
 
     ns.init(ucell, cutoff, 0);
@@ -32,32 +33,23 @@ TEST(NeighborSearchTest, TwoAtomsNeighbor)
 
     auto &list = ns.get_neighbor_list();
 
-    ASSERT_EQ(list.numneigh.size(), 2);
+    ASSERT_EQ(list.get_nlocal(), 2);
 
-    EXPECT_EQ(list.numneigh[0], 8);
-    EXPECT_EQ(list.numneigh[1], 8);
+    EXPECT_EQ(list.get_numneigh(0), 8);
+    EXPECT_EQ(list.get_numneigh(1), 8);
 }
 
 TEST(NeighborSearchTest, NoNeighbor)
 {
-    UnitCellPlus ucell;
+    ModuleBase::Matrix3 latvec;
+    latvec.e11 = 1; latvec.e12 = 0; latvec.e13 = 0;
+    latvec.e21 = 0; latvec.e22 = 1; latvec.e23 = 0;
+    latvec.e31 = 0; latvec.e32 = 0; latvec.e33 = 1;
 
-    ucell.lat0 = 1.0;
-    ucell.omega = 1.0;
-
-    ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0;
-    ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0;
-    ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1;
-
-    ucell.ntype = 1;
-    ucell.na = {2};
-    ucell.nat = 2;
-
-    ucell.tau = {
-        {0.0, 0.0, 0.0},
-        {5.0, 0.0, 0.0}
-    };
-    ucell.compute_naa();
+    UnitCellLite ucell = make_test_ucell(
+        1.0, 1.0, latvec, 1, {2},
+        {{0.0, 0.0, 0.0}, {5.0, 0.0, 0.0}}
+    );
 
     NeighborSearch ns;
 
@@ -67,80 +59,16 @@ TEST(NeighborSearchTest, NoNeighbor)
 
     auto &list = ns.get_neighbor_list();
 
-    EXPECT_EQ(list.numneigh[0], 0);
-    EXPECT_EQ(list.numneigh[1], 0);
-}
-
-TEST(NeighborSearchUnit, UCellToInputAtoms)
-{
-    UnitCellPlus ucell;
-    ucell.lat0 = 1.0;
-    ucell.omega = 1.0;
-    ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0;
-    ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0;
-    ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1;
-
-    ucell.ntype = 1;
-    ucell.na = {2};
-    ucell.nat = 2;
-    ucell.tau = {
-        {0.0, 0.0, 0.0},
-        {0.5, 0.0, 0.0}
-    };
-    ucell.compute_naa();
-    NeighborSearch ns;
-    auto inputs = ns.ucell_to_input_atoms(ucell);
-
-    EXPECT_EQ(inputs.n_atoms, 2);
-    ASSERT_EQ(inputs.InputAtom.size(), 2);
-    EXPECT_DOUBLE_EQ(inputs.InputAtom[0].position_x, 0.0);
-    EXPECT_DOUBLE_EQ(inputs.InputAtom[1].position_x, 0.5);
-    EXPECT_DOUBLE_EQ(inputs.x_low, 0.0);
-    EXPECT_DOUBLE_EQ(inputs.x_high, 0.5);
-}
-
-TEST(NeighborSearchUnit, CheckExpandAndSetMembers)
-{
-    UnitCellPlus ucell;
-    ucell.lat0 = 1.0;
-    ucell.omega = 1.0;
-    ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0;
-    ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0;
-    ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1;
-
-    ucell.ntype = 1;
-    ucell.na = {2};
-    ucell.nat = 2;
-    ucell.tau = {
-        {0.0, 0.0, 0.0},
-        {0.5, 0.0, 0.0}
-    };
-    ucell.compute_naa();
-    NeighborSearch ns;
-    ns.search_radius = 1.0; // use search radius = 1 for Check_Expand_Condition
-    ns.Check_Expand_Condition(ucell);
-
-    // For identity lattice with search_radius=1 expected ceil produce values
-    EXPECT_EQ(ns.glayerX, 2);
-    EXPECT_EQ(ns.glayerY, 2);
-    EXPECT_EQ(ns.glayerZ, 2);
-    EXPECT_EQ(ns.glayerX_minus, 1);
-
-    // Now populate all_atoms and check count
-    ns.setMemberVariables(ucell);
-    int images_x = ns.glayerX + ns.glayerX_minus; // iterations in x
-    int images_y = ns.glayerY + ns.glayerY_minus;
-    int images_z = ns.glayerZ + ns.glayerZ_minus;
-    int expected = images_x * images_y * images_z * 2; // 2 atoms per cell
-    EXPECT_EQ(static_cast<int>(ns.all_atoms.size()), expected);
+    EXPECT_EQ(list.get_numneigh(0), 0);
+    EXPECT_EQ(list.get_numneigh(1), 0);
 }
 
 TEST(NeighborSearchUnit, DistanceBox)
 {
     NeighborSearch ns;
     // set a single cell region at x=0..1,y=0..1,z=0..1
-    ns.x = 0; ns.y = 0; ns.z = 0;
-    ns.wide_x = 1.0; ns.wide_y = 1.0; ns.wide_z = 1.0;
+    ns.set_position(0, 0, 0);
+    ns.set_width(1.0, 1.0, 1.0);
 
     double inside = ns.distance(0.2, 0.5, 0.5, 0.0, 0.0, 0.0);
     EXPECT_DOUBLE_EQ(inside, 0.0);
@@ -169,32 +97,6 @@ TEST(NeighborSearchUnit, DecomposeCases)
     EXPECT_EQ(nz, 7);
 }
 
-TEST(NeighborSearchUnit, UCellToInputAtomsMultipleTypes)
-{
-    UnitCellPlus ucell;
-    ucell.lat0 = 1.0;
-    ucell.omega = 1.0;
-    ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0;
-    ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0;
-    ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1;
-
-    ucell.ntype = 2;
-    ucell.na = {1, 2};
-    ucell.nat = 3;
-    ucell.tau = {
-        {0.0, 0.0, 0.0},
-        {0.5, 0.0, 0.0},
-        {0.0, 0.5, 0.0}
-    };
-    ucell.compute_naa();
-    NeighborSearch ns;
-    auto inputs = ns.ucell_to_input_atoms(ucell);
-
-    EXPECT_EQ(inputs.n_atoms, 3);
-    ASSERT_EQ(inputs.InputAtom.size(), 3);
-    EXPECT_DOUBLE_EQ(inputs.InputAtom[2].position_y, 0.5);
-}
-
 TEST(NeighborSearchUnit, DecomposePrimeNumber)
 {
     NeighborSearch ns;
@@ -208,83 +110,70 @@ TEST(NeighborSearchUnit, DecomposePrimeNumber)
 
 TEST(NeighborSearchUnit, NonOrthogonalLatticeExpand)
 {
-    UnitCellPlus ucell;
-    ucell.lat0 = 1.0;
-    ucell.omega = 1.0;
+    ModuleBase::Matrix3 latvec;
     // skewed lattice
-    ucell.latvec.e11 = 1; ucell.latvec.e12 = 0.3; ucell.latvec.e13 = 0.0;
-    ucell.latvec.e21 = 0.1; ucell.latvec.e22 = 1.0; ucell.latvec.e23 = 0.0;
-    ucell.latvec.e31 = 0.0; ucell.latvec.e32 = 0.0; ucell.latvec.e33 = 1.0;
+    latvec.e11 = 1; latvec.e12 = 0.3; latvec.e13 = 0.0;
+    latvec.e21 = 0.1; latvec.e22 = 1.0; latvec.e23 = 0.0;
+    latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = 1.0;
 
-    ucell.ntype = 1;
-    ucell.na = {1};
-    ucell.nat = 1;
-    ucell.tau = {{0.0, 0.0, 0.0}};
-    ucell.compute_naa();
+    UnitCellLite ucell = make_test_ucell(
+        1.0, 1.0, latvec, 1, {1},
+        {{0.0, 0.0, 0.0}}
+    );
 
     NeighborSearch ns;
-    ns.search_radius = 2.5;
-    ns.Check_Expand_Condition(ucell);
+    ns.init(ucell, 2.5, 0);
     // for skewed lattice, expansion layers should be >= 1
-    EXPECT_GE(ns.glayerX, 1);
-    EXPECT_GE(ns.glayerY, 1);
-    EXPECT_GE(ns.glayerZ, 1);
+    EXPECT_GE(ns.get_glayerX(), 1);
+    EXPECT_GE(ns.get_glayerY(), 1);
+    EXPECT_GE(ns.get_glayerZ(), 1);
 }
 
-// --- additional tests to cover remaining branches in neighbor_search.cpp ---
-
 TEST(NeighborSearchInit_WideZero_CentralInside, SingleAtomCell)
 {
-    UnitCellPlus ucell;
-    ucell.lat0 = 1.0;
-    ucell.omega = 1.0;
-    ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0;
-    ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0;
-    ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1;
-
-    ucell.ntype = 1;
-    ucell.na = {1};
-    ucell.nat = 1;
-    ucell.tau = {{0.0, 0.0, 0.0}};
-    ucell.compute_naa();
+    ModuleBase::Matrix3 latvec;
+    latvec.e11 = 1; latvec.e12 = 0; latvec.e13 = 0;
+    latvec.e21 = 0; latvec.e22 = 1; latvec.e23 = 0;
+    latvec.e31 = 0; latvec.e32 = 0; latvec.e33 = 1;
+
+    UnitCellLite ucell = make_test_ucell(
+        1.0, 1.0, latvec, 1, {1},
+        {{0.0, 0.0, 0.0}}
+    );
 
     NeighborSearch ns;
     // choose sr small enough; with mpi_size fixed to 1 in init, wide_* become 0
     ns.init(ucell, 0.1, 0);
     // central cell atom should be counted as inside
-    EXPECT_EQ(ns.inside_atoms.size(), 1);
-    EXPECT_EQ(ns.neighbor_list.nlocal, static_cast<int>(ns.inside_atoms.size()));
+    EXPECT_EQ(ns.get_inside_atoms().size(), 1);
+    EXPECT_EQ(ns.get_neighbor_list().get_nlocal(), static_cast<int>(ns.get_inside_atoms().size()));
 }
 
 TEST(NeighborSearchInit_MpiRankIndexing, RankValues)
 {
-    UnitCellPlus ucell;
-    ucell.lat0 = 1.0;
-    ucell.omega = 1.0;
-    ucell.latvec.e11 = 1; ucell.latvec.e12 = 0; ucell.latvec.e13 = 0;
-    ucell.latvec.e21 = 0; ucell.latvec.e22 = 1; ucell.latvec.e23 = 0;
-    ucell.latvec.e31 = 0; ucell.latvec.e32 = 0; ucell.latvec.e33 = 1;
-
-    ucell.ntype = 1;
-    ucell.na = {1};
-    ucell.nat = 1;
-    ucell.tau = {{0.0, 0.0, 0.0}};
-    ucell.compute_naa();
+    ModuleBase::Matrix3 latvec;
+    latvec.e11 = 1; latvec.e12 = 0; latvec.e13 = 0;
+    latvec.e21 = 0; latvec.e22 = 1; latvec.e23 = 0;
+    latvec.e31 = 0; latvec.e32 = 0; latvec.e33 = 1;
+
+    UnitCellLite ucell = make_test_ucell(
+        1.0, 1.0, latvec, 1, {1},
+        {{0.0, 0.0, 0.0}}
+    );
 
     NeighborSearch ns0;
     ns0.init(ucell, 0.5, 0);
     // with mpi_size fixed to 1 in init, nx=ny=nz=1; for rank 0 expect x=y=0,z=0
-    EXPECT_EQ(ns0.x, 0);
-    EXPECT_EQ(ns0.y, 0);
-    EXPECT_EQ(ns0.z, 0);
-
+    EXPECT_EQ(ns0.get_x(), 0);
+    EXPECT_EQ(ns0.get_y(), 0);
+    EXPECT_EQ(ns0.get_z(), 0);
 }
 
 TEST(NeighborSearchDistance_OutsideCases, VariousAxes)
 {
     NeighborSearch ns;
-    ns.x = 0; ns.y = 0; ns.z = 0;
-    ns.wide_x = 2.0; ns.wide_y = 3.0; ns.wide_z = 4.0;
+    ns.set_position(0, 0, 0);
+    ns.set_width(2.0, 3.0, 4.0);
 
     // position inside box along x (no dx), but outside along y by above high bound
     double d = ns.distance(0.5, 4.5, 1.0, 0.0, 0.0, 0.0);
@@ -313,4 +202,33 @@ TEST(NeighborSearchDecompose_SmallSizes, TwoAndOne)
     EXPECT_EQ(nz, 1);
 }
 
-// end of additional tests
+TEST(NeighborSearchUnit, ExpansionLayersAndAtomCount)
+{
+    ModuleBase::Matrix3 latvec;
+    latvec.e11 = 1; latvec.e12 = 0; latvec.e13 = 0;
+    latvec.e21 = 0; latvec.e22 = 1; latvec.e23 = 0;
+    latvec.e31 = 0; latvec.e32 = 0; latvec.e33 = 1;
+
+    UnitCellLite ucell = make_test_ucell(
+        1.0, 1.0, latvec, 1, {2},
+        {{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}
+    );
+
+    NeighborSearch ns;
+    ns.init(ucell, 1.0, 0);
+
+    // For identity lattice with search_radius=1 expected ceil produce values
+    EXPECT_EQ(ns.get_glayerX(), 2);
+    EXPECT_EQ(ns.get_glayerY(), 2);
+    EXPECT_EQ(ns.get_glayerZ(), 2);
+    EXPECT_EQ(ns.get_glayerX_minus(), 1);
+
+    // Check atom count
+    int images_x = ns.get_glayerX() + ns.get_glayerX_minus();
+    int images_y = ns.get_glayerY() + ns.get_glayerY_minus();
+    int images_z = ns.get_glayerZ() + ns.get_glayerZ_minus();
+    int expected = images_x * images_y * images_z * 2; // 2 atoms per cell
+    EXPECT_EQ(static_cast<int>(ns.get_all_atoms().size()), expected);
+}
+
+// end of additional tests
\ No newline at end of file
diff --git a/source/source_cell/module_neighlist/unitcell_interface.h b/source/source_cell/module_neighlist/unitcell_interface.h
deleted file mode 100644
index 25610d471b3..00000000000
--- a/source/source_cell/module_neighlist/unitcell_interface.h
+++ /dev/null
@@ -1,25 +0,0 @@
-#ifndef UNITCELL_INTERFACE_H
-#define UNITCELL_INTERFACE_H
-
-
-#include "source_base/vector3.h"
-#include "source_base/matrix3.h"
-
-
-class IAtomProvider {
-public:
-    virtual ~IAtomProvider() = default;
-
-    virtual double get_lat0() const = 0;
-    virtual double get_omega() const = 0;
-    virtual const ModuleBase::Matrix3& get_latvec() const = 0;
-
-
- 
-    virtual int get_natom() const = 0;
-    virtual int get_na(int i) const = 0;
-    virtual int get_ntype() const = 0;
-    virtual ModuleBase::Vector3<double> get_tauu(int i,int j) const = 0;
-};
-
-#endif // UNITCELL_INTERFACE_H
\ No newline at end of file
diff --git a/source/source_cell/module_neighlist/unitcell_lite.cpp b/source/source_cell/module_neighlist/unitcell_lite.cpp
new file mode 100644
index 00000000000..8475f81cf29
--- /dev/null
+++ b/source/source_cell/module_neighlist/unitcell_lite.cpp
@@ -0,0 +1,92 @@
+#include "unitcell_lite.h"
+
+#include <cassert>
+
+// === AtomProvider interface implementation ===
+
+double UnitCellLite::get_lat0() const {
+    return lat0_;
+}
+
+double UnitCellLite::get_omega() const {
+    return omega_;
+}
+
+const ModuleBase::Matrix3& UnitCellLite::get_latvec() const {
+    return latvec_;
+}
+
+int UnitCellLite::get_natom() const {
+    return nat_;
+}
+
+int UnitCellLite::get_na(int i) const {
+    assert(i >= 0 && i < ntype_);
+    return na_[i];
+}
+
+int UnitCellLite::get_ntype() const {
+    return ntype_;
+}
+
+ModuleBase::Vector3<double> UnitCellLite::get_tau(int i, int j) const {
+    assert(i >= 0 && i < ntype_);
+    assert(j >= 0 && j < na_[i]);
+    if (i == 0) {
+        return tau_[j];
+    }
+    return tau_[naa_[i - 1] + j];
+}
+
+// === Setter methods ===
+
+void UnitCellLite::set_lat0(double lat0) {
+    lat0_ = lat0;
+}
+
+void UnitCellLite::set_omega(double omega) {
+    omega_ = omega;
+}
+
+void UnitCellLite::set_latvec(const ModuleBase::Matrix3& latvec) {
+    latvec_ = latvec;
+}
+
+void UnitCellLite::set_lattice(double lat0, double omega, const ModuleBase::Matrix3& latvec) {
+    lat0_ = lat0;
+    omega_ = omega;
+    latvec_ = latvec;
+}
+
+void UnitCellLite::set_atoms(int ntype,
+                              const std::vector<int>& na,
+                              const std::vector<ModuleBase::Vector3<double>>& tau) {
+    assert(ntype >= 0);
+    assert(na.size() == static_cast<size_t>(ntype));
+    
+    ntype_ = ntype;
+    na_ = na;
+    tau_ = tau;
+    
+    // compute total number of atoms
+    nat_ = 0;
+    for (int i = 0; i < ntype_; ++i) {
+        nat_ += na_[i];
+    }
+    assert(tau_.size() == static_cast<size_t>(nat_));
+    
+    // compute cumulative counts
+    compute_naa_();
+}
+
+// === Internal methods ===
+
+void UnitCellLite::compute_naa_() {
+    naa_.resize(na_.size());
+    if (naa_.size() > 0) {
+        naa_[0] = na_[0];
+    }
+    for (size_t i = 1; i < naa_.size(); ++i) {
+        naa_[i] = naa_[i - 1] + na_[i];
+    }
+}
\ No newline at end of file
diff --git a/source/source_cell/module_neighlist/unitcell_lite.h b/source/source_cell/module_neighlist/unitcell_lite.h
new file mode 100644
index 00000000000..e951945799f
--- /dev/null
+++ b/source/source_cell/module_neighlist/unitcell_lite.h
@@ -0,0 +1,169 @@
+#ifndef UNITCELL_LITE_H
+#define UNITCELL_LITE_H
+
+#include "source_cell/module_neighlist/atom_provider.h"
+#include <vector>
+
+/**
+ * @brief A lightweight unit cell class for molecular dynamics simulations.
+ *
+ * This class provides a minimal set of unit cell information needed for
+ * large-scale molecular dynamics simulations (e.g., billion-atom simulations).
+ * It implements the AtomProvider interface and stores only essential data:
+ * lattice parameters and atomic coordinates.
+ *
+ * Compared to the full UnitCell class, UnitCellLite has significantly lower
+ * memory overhead by omitting electronic structure-related data such as
+ * pseudopotentials, orbitals, magnetism, and symmetry information.
+ *
+ * @see AtomProvider
+ * @see UnitCell
+ */
+class UnitCellLite : public AtomProvider
+{
+public:
+    /**
+     * @brief Default constructor.
+     *
+     * Initializes all data members to zero/empty state.
+     */
+    UnitCellLite() = default;
+
+    /**
+     * @brief Default destructor.
+     */
+    ~UnitCellLite() = default;
+
+    // ========== AtomProvider interface implementation ==========
+
+    /**
+     * @brief Get the lattice constant in Bohr.
+     * @return Lattice constant lat0.
+     */
+    double get_lat0() const override;
+
+    /**
+     * @brief Get the unit cell volume.
+     * @return Cell volume omega in Bohr^3.
+     */
+    double get_omega() const override;
+
+    /**
+     * @brief Get the lattice vectors.
+     * @return Reference to the 3x3 matrix of lattice vectors.
+     */
+    const ModuleBase::Matrix3& get_latvec() const override;
+
+    /**
+     * @brief Get the total number of atoms.
+     * @return Total atom count nat.
+     */
+    int get_natom() const override;
+
+    /**
+     * @brief Get the number of atoms for a given type.
+     * @param i Atom type index (0-based).
+     * @return Number of atoms of type i.
+     * @note Asserts that i is in valid range [0, ntype_).
+     */
+    int get_na(int i) const override;
+
+    /**
+     * @brief Get the number of atom types.
+     * @return Number of atom types ntype.
+     */
+    int get_ntype() const override;
+
+    /**
+     * @brief Get the coordinate of atom (type i, index j).
+     * @param i Atom type index (0-based).
+     * @param j Atom index within type i (0-based).
+     * @return Cartesian coordinate of the atom in Bohr.
+     * @note Asserts that i and j are in valid ranges.
+     */
+    ModuleBase::Vector3<double> get_tau(int i, int j) const override;
+
+    // ========== Setter methods ==========
+
+    /**
+     * @brief Set the lattice constant.
+     * @param lat0 Lattice constant in Bohr.
+     */
+    void set_lat0(double lat0);
+
+    /**
+     * @brief Set the unit cell volume.
+     * @param omega Cell volume in Bohr^3.
+     */
+    void set_omega(double omega);
+
+    /**
+     * @brief Set the lattice vectors.
+     * @param latvec 3x3 matrix of lattice vectors.
+     */
+    void set_latvec(const ModuleBase::Matrix3& latvec);
+
+    /**
+     * @brief Set all lattice parameters together.
+     * @param lat0 Lattice constant in Bohr.
+     * @param omega Cell volume in Bohr^3.
+     * @param latvec 3x3 matrix of lattice vectors.
+     */
+    void set_lattice(double lat0, double omega, const ModuleBase::Matrix3& latvec);
+
+    /**
+     * @brief Set atom information for all types.
+     *
+     * This method sets the number of atom types, the count of atoms per type,
+     * and all atomic coordinates. It automatically computes the total atom
+     * count (nat_) and the cumulative atom counts (naa_).
+     *
+     * @param ntype Number of atom types.
+     * @param na Vector of atom counts for each type [ntype].
+     * @param tau Vector of all atomic coordinates [nat].
+     *
+     * @note Asserts that na.size() == ntype and tau.size() == sum(na).
+     */
+    void set_atoms(int ntype,
+                   const std::vector<int>& na,
+                   const std::vector<ModuleBase::Vector3<double>>& tau);
+
+private:
+    // ========== Data members ==========
+
+    /// Lattice constant in Bohr
+    double lat0_ = 0.0;
+
+    /// Unit cell volume in Bohr^3
+    double omega_ = 0.0;
+
+    /// Total number of atoms
+    int nat_ = 0;
+
+    /// Number of atom types
+    int ntype_ = 0;
+
+    /// Lattice vectors (3x3 matrix)
+    ModuleBase::Matrix3 latvec_;
+
+    /// Number of atoms for each type [ntype]
+    std::vector<int> na_;
+
+    /// Cumulative sum of na: naa_[i] = na_[0] + na_[1] + ... + na_[i]
+    std::vector<int> naa_;
+
+    /// Atomic coordinates in Cartesian (Bohr) [nat]
+    std::vector<ModuleBase::Vector3<double>> tau_;
+
+    // ========== Internal methods ==========
+
+    /**
+     * @brief Compute cumulative atom counts from na_.
+     *
+     * Updates naa_ such that naa_[i] = sum of na_[0] to na_[i].
+     * Called internally by set_atoms().
+     */
+    void compute_naa_();
+};
+
+#endif // UNITCELL_LITE_H
\ No newline at end of file
diff --git a/source/source_cell/module_neighlist/unitcell_plus.h b/source/source_cell/module_neighlist/unitcell_plus.h
deleted file mode 100644
index 3e15fed210d..00000000000
--- a/source/source_cell/module_neighlist/unitcell_plus.h
+++ /dev/null
@@ -1,73 +0,0 @@
-#pragma once
-
-#include "source_cell/module_neighlist/unitcell_interface.h"
-#include <vector>
-
-class UnitCellPlus : public IAtomProvider
-{
-public:
-    UnitCellPlus()=default;
-    ~UnitCellPlus()=default;
-
-    double get_lat0() const override {
-        return lat0;
-    }
-
-    double get_omega() const override {
-        return omega;
-    }
-
-    const ModuleBase::Matrix3& get_latvec() const override {
-        return latvec;
-    }
-
-    int get_natom() const override {
-        return nat;
-    }
-
-    int get_na(int i) const override {
-        assert(i >= 0 && i < ntype);
-        return na[i];
-    }
-
-    int get_ntype() const override {
-        return ntype;
-    }
-
-    ModuleBase::Vector3<double> get_tauu(int i, int j) const override {
-        assert(i >= 0 && i < ntype);
-        assert(j >= 0 && j < na[i]);
-        //assert(naa.size() > 0);
-        if(i==0)
-        {
-            return tau[j];
-        }
-        return tau[naa[i-1]+j];
-    }
-
-    double lat0;
-    double omega;
-    int nat;
-    std::vector<int> na;
-    // naa is the cumulative sum of na: naa[i] = na[0] + na[1] + ... + na[i]
-    std::vector<int> naa;
-    int ntype;
-    ModuleBase::Matrix3 latvec;
-    std::vector<ModuleBase::Vector3<double>> tau;
-
-    // Compute cumulative counts in naa from current na vector.
-    // Assumption: naa[i] == sum_{k=0..i} na[k]. Call this after na is set/modified.
-    void compute_naa()
-    {
-        naa.resize(na.size());
-        if(naa.size()>0)
-        {
-            naa[0] = na[0];
-        }
-        for(int i=1;i<naa.size();i++)
-        {
-            naa[i] = naa[i-1] + na[i];
-        }
-    }
-
-};
\ No newline at end of file
diff --git a/source/source_cell/unitcell.h b/source/source_cell/unitcell.h
index 0b62d0ae5b3..df426712b47 100644
--- a/source/source_cell/unitcell.h
+++ b/source/source_cell/unitcell.h
@@ -5,14 +5,14 @@
 #include "source_cell/sep_cell.h"
 #include "source_estate/magnetism.h"
 #include "module_symmetry/symmetry.h"
-#include "source_cell/module_neighlist/unitcell_interface.h"
+#include "source_cell/module_neighlist/atom_provider.h"
 
 #ifdef __LCAO
 #include "setup_nonlocal.h"
 #endif
 
 // provide the basic information about unitcell.
-class UnitCell : public IAtomProvider {
+class UnitCell : public AtomProvider {
   public:
     double get_lat0() const override {
         return lat0;
@@ -38,7 +38,7 @@ class UnitCell : public IAtomProvider {
         return ntype;
     }
 
-    ModuleBase::Vector3<double> get_tauu(int i, int j) const override {
+    ModuleBase::Vector3<double> get_tau(int i, int j) const override {
         return atoms[i].tau[j];
     }
 
diff --git a/source/source_esolver/esolver_lj.cpp b/source/source_esolver/esolver_lj.cpp
index 72db00a1db2..3ddba86adc2 100644
--- a/source/source_esolver/esolver_lj.cpp
+++ b/source/source_esolver/esolver_lj.cpp
@@ -11,27 +11,25 @@
 namespace ModuleESolver
 {
 
-    UnitCellPlus ESolver_LJ::change_from_ucell_to_ucell_plus(const UnitCell& ucell)
+    UnitCellLite ESolver_LJ::change_from_ucell_to_ucell_lite(const UnitCell& ucell)
     {
-        UnitCellPlus ucell_plus;
-        ucell_plus.lat0 = ucell.lat0;
-        ucell_plus.omega = ucell.omega;
-        ucell_plus.nat = ucell.nat;
-        for(int i=0;i<ucell.ntype;i++)
-        {
-            ucell_plus.na.push_back(ucell.atoms[i].na);
-        }
-        ucell_plus.ntype = ucell.ntype;
-        ucell_plus.latvec = ucell.latvec;
-        for(int i=0;i<ucell.ntype;i++)
-        {
-            for(int j=0;j<ucell.atoms[i].na;j++)
-            {
-                ucell_plus.tau.push_back(ucell.atoms[i].tau[j]);
+        UnitCellLite ucell_lite;
+        
+        // Set lattice parameters
+        ucell_lite.set_lattice(ucell.lat0, ucell.omega, ucell.latvec);
+        
+        // Build atom information
+        std::vector<int> na;
+        std::vector<ModuleBase::Vector3<double>> tau;
+        for (int i = 0; i < ucell.ntype; i++) {
+            na.push_back(ucell.atoms[i].na);
+            for (int j = 0; j < ucell.atoms[i].na; j++) {
+                tau.push_back(ucell.atoms[i].tau[j]);
             }
         }
-        ucell_plus.compute_naa();
-        return ucell_plus;
+        ucell_lite.set_atoms(ucell.ntype, na, tau);
+        
+        return ucell_lite;
     }
 
 void ESolver_LJ::before_all_runners(UnitCell& ucell, const Input_para& inp)
@@ -57,9 +55,9 @@ void ESolver_LJ::before_all_runners(UnitCell& ucell, const Input_para& inp)
 
 void ESolver_LJ::runner(UnitCell& ucell, const int istep)
 {
-    UnitCellPlus ucell_plus = change_from_ucell_to_ucell_plus(ucell);
+    UnitCellLite ucell_lite = change_from_ucell_to_ucell_lite(ucell);
     NeighborSearch neighbor_search;
-    neighbor_search.init(ucell_plus, search_radius, 0);
+    neighbor_search.init(ucell_lite, search_radius, 0);
     neighbor_search.build_neighbors();
 
     double distance = 0.0;
@@ -71,18 +69,21 @@ void ESolver_LJ::runner(UnitCell& ucell, const int istep)
     lj_virial.zero_out();
 
     ModuleBase::Vector3<double> tau1, tau2, dtau;
+    const NeighborList& neighbor_list = neighbor_search.get_neighbor_list();
+    const std::vector<NeighborAtom>& all_atoms = neighbor_search.get_all_atoms();
     for (int it = 0; it < ucell.ntype; ++it)
     {
         Atom* atom1 = &ucell.atoms[it];
         for (int ia = 0; ia < atom1->na; ++ia)
         {
             tau1 = atom1->tau[ia];
-            for (int ad = 0; ad < neighbor_search.neighbor_list.numneigh[index]; ++ad)
+            for (int ad = 0; ad < neighbor_list.get_numneigh(index); ++ad)
             {
-                tau2.x = neighbor_search.all_atoms[neighbor_search.neighbor_list.firstneigh[index][ad]].position_x;
-                tau2.y = neighbor_search.all_atoms[neighbor_search.neighbor_list.firstneigh[index][ad]].position_y;
-                tau2.z = neighbor_search.all_atoms[neighbor_search.neighbor_list.firstneigh[index][ad]].position_z;
-                int it2 = neighbor_search.all_atoms[neighbor_search.neighbor_list.firstneigh[index][ad]].atom_type;
+                const NeighborAtom& neighbor_atom = all_atoms[neighbor_list.get_firstneigh(index)[ad]];
+                tau2.x = neighbor_atom.position_x;
+                tau2.y = neighbor_atom.position_y;
+                tau2.z = neighbor_atom.position_z;
+                int it2 = neighbor_atom.atom_type;
                 dtau = (tau1 - tau2) * ucell.lat0;
                 distance = dtau.norm();
                 if (distance < lj_rcut(it, it2))
diff --git a/source/source_esolver/esolver_lj.h b/source/source_esolver/esolver_lj.h
index a8d8d9d1787..1a96510eaa4 100644
--- a/source/source_esolver/esolver_lj.h
+++ b/source/source_esolver/esolver_lj.h
@@ -2,7 +2,7 @@
 #define ESOLVER_LJ_H
 
 #include "esolver.h"
-#include "source_cell/module_neighlist/unitcell_plus.h"
+#include "source_cell/module_neighlist/unitcell_lite.h"
 
 namespace ModuleESolver
 {
@@ -15,7 +15,7 @@ namespace ModuleESolver
             classname = "ESolver_LJ";
         }
 
-        UnitCellPlus change_from_ucell_to_ucell_plus(const UnitCell& ucell);
+        UnitCellLite change_from_ucell_to_ucell_lite(const UnitCell& ucell);
 
         void before_all_runners(UnitCell& ucell, const Input_para& inp) override;
 
diff --git a/source/source_md/test/CMakeLists.txt b/source/source_md/test/CMakeLists.txt
index c4e3a1c2d2f..3ff705c4f4a 100644
--- a/source/source_md/test/CMakeLists.txt
+++ b/source/source_md/test/CMakeLists.txt
@@ -47,6 +47,8 @@ list(APPEND depend_files
   ../../source_cell/module_neighbor/sltk_grid_driver.cpp
   ../../source_cell/module_neighlist/neighbor_search.cpp
   ../../source_cell/module_neighlist/bin_manager.cpp
+  ../../source_cell/module_neighlist/page_allocator.cpp
+  ../../source_cell/module_neighlist/unitcell_lite.cpp
   ../../source_io/module_output/output.cpp
   ../../source_io/module_output/output_log.cpp
   ../../source_io/module_output/print_info.cpp