diff --git a/source/source_base/kernels/cuda/sph_harm_gpu.cuh b/source/source_base/kernels/cuda/sph_harm_gpu.cuh
index d4fa5f5666..caf29795d6 100644
--- a/source/source_base/kernels/cuda/sph_harm_gpu.cuh
+++ b/source/source_base/kernels/cuda/sph_harm_gpu.cuh
@@ -4,36 +4,15 @@
 
 namespace ModuleBase {
 
-/// Spherical harmonics computation (table lookup method)
-/// Directly uses constexpr ylmcoef, compiler auto-inlines
-/// @param nwl Maximum angular momentum L (0 <= nwl <= 5)
-/// @param x,y,z Direction vector (need not be normalized, normalization is done internally)
-/// @param ylma Output array, size (nwl+1)^2
-__device__ static void sph_harm(
+/// Evaluate the existing spherical-harmonic recurrence directly.
+/// This helper performs no input normalization and no zero-vector fallback.
+__device__ static void sph_harm_direct(
     const int nwl,
-    const double x_in,
-    const double y_in,
-    const double z_in,
+    const double x,
+    const double y,
+    const double z,
     double* __restrict__ ylma)
 {
-   // Normalize the input direction vector
-   double r = sqrt(x_in * x_in + y_in * y_in + z_in * z_in);
-   double x, y, z;
-   if (r < 1e-10)
-   {
-       // At origin, default to z-axis direction
-       x = 0.0;
-       y = 0.0;
-       z = 1.0;
-   }
-   else
-   {
-       const double inv_r = 1.0 / r;
-       x = x_in * inv_r;
-       y = y_in * inv_r;
-       z = z_in * inv_r;
-   }
-
    /***************************
    L = 0
    ***************************/
@@ -147,6 +126,39 @@ __device__ static void sph_harm(
        return;
 }
 
+/// Spherical harmonics computation (table lookup method)
+/// Directly uses constexpr ylmcoef, compiler auto-inlines
+/// @param nwl Maximum angular momentum L (0 <= nwl <= 5)
+/// @param x,y,z Direction vector (need not be normalized, normalization is done internally)
+/// @param ylma Output array, size (nwl+1)^2
+__device__ static void sph_harm(
+    const int nwl,
+    const double x_in,
+    const double y_in,
+    const double z_in,
+    double* __restrict__ ylma)
+{
+   // Normalize the input direction vector
+   double r = sqrt(x_in * x_in + y_in * y_in + z_in * z_in);
+   double x, y, z;
+   if (r < 1e-10)
+   {
+       // At origin, default to z-axis direction
+       x = 0.0;
+       y = 0.0;
+       z = 1.0;
+   }
+   else
+   {
+       const double inv_r = 1.0 / r;
+       x = x_in * inv_r;
+       y = y_in * inv_r;
+       z = z_in * inv_r;
+   }
+
+    sph_harm_direct(nwl, x, y, z, ylma);
+}
+
 /// Spherical harmonics and gradient computation
 __device__ static void grad_rl_sph_harm(
     const int nwl,
diff --git a/source/source_hamilt/module_gint/kernel/phi_operator_kernel.cuh b/source/source_hamilt/module_gint/kernel/phi_operator_kernel.cuh
index 7bc7c70594..00a53860c7 100644
--- a/source/source_hamilt/module_gint/kernel/phi_operator_kernel.cuh
+++ b/source/source_hamilt/module_gint/kernel/phi_operator_kernel.cuh
@@ -53,6 +53,11 @@ __global__ void set_phi_kernel(
         const double3 coord = make_double3(mgrid_pos.x-rcoord.x,                    // coord is the relative coordinate of an atom and a meshgrid
                                            mgrid_pos.y-rcoord.y,
                                            mgrid_pos.z-rcoord.z);
+        // Preserve the existing near-origin behavior. Only the exact
+        // atomic grid point follows the CPU direct-recurrence semantics.
+        const bool exact_origin
+            = (coord.x == 0.0 && coord.y == 0.0 && coord.z == 0.0);
+
         double dist = norm3d(coord.x, coord.y, coord.z);
         if (dist < rcut[atom_type])
         {
@@ -61,7 +66,20 @@ __global__ void set_phi_kernel(
             // since nwl is less or equal than 5, the size of ylma is (5+1)^2
             double ylma[36];
             const int nwl = ucell_atom_nwl[atom_type];
-            sph_harm(nwl, coord.x/dist, coord.y/dist, coord.z/dist, ylma);
+            if (exact_origin)
+            {
+                ModuleBase::sph_harm_direct(
+                    nwl, 0.0, 0.0, 0.0, ylma);
+            }
+            else
+            {
+                sph_harm(
+                    nwl,
+                    coord.x/dist,
+                    coord.y/dist,
+                    coord.z/dist,
+                    ylma);
+            }
 
             const double pos = dist / dr_uniform;
             const int ip = static_cast<int>(pos);
diff --git a/tests/03_NAO_multik/CASES_GPU.txt b/tests/03_NAO_multik/CASES_GPU.txt
index 88d9f4478a..1453c8cf28 100644
--- a/tests/03_NAO_multik/CASES_GPU.txt
+++ b/tests/03_NAO_multik/CASES_GPU.txt
@@ -33,7 +33,7 @@ scf_out_hsk
 scf_out_hsk_binary
 scf_out_hsr
 scf_out_hsr_binary_spin2
-#scf_out_hsr_spin4
+scf_out_hsr_spin4
 scf_out_dh_t
 scf_out_dos_spin4
 scf_out_mul
@@ -46,7 +46,7 @@ nscf_out_dos
 nscf_out_band_pband
 nscf_out_pot1
 nscf_out_mul
-#nscf_out_hsr_tr_rr
+nscf_out_hsr_tr_rr
 relax_bfgs2
 relax_old_cg
 relax_cell