diff --git a/docs/scanl_laplacian_implementation_report.md b/docs/scanl_laplacian_implementation_report.md
new file mode 100644
index 00000000000..35b5d493ffc
--- /dev/null
+++ b/docs/scanl_laplacian_implementation_report.md
@@ -0,0 +1,333 @@
+# SCANL 有限差分 Laplacian 实现
+
+## 1. 问题
+
+SCANL (deorbitalized SCAN) 将 SCAN 对动能密度 τ 的依赖替换为对密度 Laplacian ∇²ρ 的依赖。其 Kohn-Sham 势包含项
+
+$$V_{xc}^{\nabla^2\rho}(\mathbf{r}) = e^2\,\nabla^2\!\left(\frac{\partial(\rho\varepsilon_{xc})}{\partial(\nabla^2\rho)}\right) \equiv e^2\,\nabla^2(\text{vlapl})$$
+
+在平面波基组下，倒易空间中 ∇² → −|G|²，因此
+
+$$\tilde{V}_{xc}^{\nabla^2\rho}(\mathbf{G}) = -e^2\,|\mathbf{G}|^2\,\widetilde{\text{vlapl}}(\mathbf{G})$$
+
+|G|² 无界增长导致 SCF 发散：vlapl(G) 的高频数值噪声被 |G|² 放大后反馈到密度，形成正反馈。此发散不依赖 mixing 策略，是 Laplace 逆变换的数学病态性。
+
+---
+
+## 2. 有限差分 Laplacian
+
+### 2.1 核心思想
+
+FFT 网格上的密度是离散函数。对离散函数求 Laplacian，有限差分 (FD) 是自然且唯一自洽的方法。FD Laplacian 的 G 空间 kernel 在低 G 时与谱 kernel 一致，在高 G 时有界衰减，不会放大高频噪声。
+
+### 2.2 坐标系与度规
+
+设晶格矢量 $\mathbf{a}_\alpha$（α = 1,2,3），以 lat0 为单位存储在 `rhopw->latvec` 中。分数坐标
+
+$$u^\alpha = n_\alpha / N_\alpha, \quad n_\alpha = 0, 1, \ldots, N_\alpha - 1$$
+
+其中 $N_\alpha$ = `rhopw->nx, ny, nz`。实空间位置
+
+$$\mathbf{r} = \text{lat0} \sum_\alpha u^\alpha \mathbf{a}_\alpha$$
+
+Laplacian 在分数坐标下为
+
+$$\nabla^2 = \sum_{\alpha,\beta} g^{\alpha\beta} \frac{\partial^2}{\partial u^\alpha \partial u^\beta}$$
+
+其中逆变度规张量
+
+$$g^{\alpha\beta} = \frac{\text{GGT}[\alpha][\beta]}{\text{lat0}^2}$$
+
+GGT = G · Gᵀ 是倒格子度规矩阵（ABACUS 中的 `rhopw->GGT`），G = (latvec)⁻ᵀ。
+
+### 2.3 平面波的谱 Laplacian
+
+对平面波 $f(\mathbf{r}) = \exp(2\pi i \sum_\gamma m_\gamma u^\gamma)$，其中 $m_\gamma$ 为 Miller 指数（`rhopw->gdirect[ig]`）：
+
+$$\frac{\partial^2 f}{\partial u^\alpha \partial u^\beta}\bigg|_{\text{spectral}} = -(2\pi)^2 m_\alpha m_\beta \, f$$
+
+组合得谱 Laplacian：
+
+$$\nabla^2_{\text{spectral}} f = -\underbrace{\left[\sum_{\alpha,\beta} \text{GGT}[\alpha][\beta] \cdot m_\alpha m_\beta\right]}_{\text{gg}_{\text{spectral}}} \cdot \text{tpiba2} \cdot f$$
+
+其中 tpiba2 = (2π/lat0)²。这就是 ABACUS 中 `gg[ig] * tpiba2` 的来源。
+
+### 2.4 平面波的 FD Laplacian
+
+在分数坐标的均匀网格上对二阶导数做中心差分：
+
+**对角项** (α = β)：
+
+$$\frac{\partial^2 f}{\partial (u^\alpha)^2}\bigg|_{\text{FD}} = \frac{f(u^\alpha + 1/N_\alpha) + f(u^\alpha - 1/N_\alpha) - 2f(u^\alpha)}{(1/N_\alpha)^2}$$
+
+对平面波：
+
+$$= -2N_\alpha^2 \left(1 - \cos\frac{2\pi m_\alpha}{N_\alpha}\right) f$$
+
+**交叉项** (α ≠ β)：
+
+$$\frac{\partial^2 f}{\partial u^\alpha \partial u^\beta}\bigg|_{\text{FD}} = \frac{1}{4 N_\alpha N_\beta}\Big[f\Big(\frac{+1}{N_\alpha}, \frac{+1}{N_\beta}\Big) - f\Big(\frac{+1}{N_\alpha}, \frac{-1}{N_\beta}\Big) - f\Big(\frac{-1}{N_\alpha}, \frac{+1}{N_\beta}\Big) + f\Big(\frac{-1}{N_\alpha}, \frac{-1}{N_\beta}\Big)\Big]$$
+
+对平面波：
+
+$$= -N_\alpha N_\beta \sin\frac{2\pi m_\alpha}{N_\alpha} \sin\frac{2\pi m_\beta}{N_\beta} \, f$$
+
+### 2.5 FD kernel
+
+定义 FD 因子矩阵：
+
+$$\text{FD}_{\alpha\alpha} = 2N_\alpha^2\left(1 - \cos\frac{2\pi m_\alpha}{N_\alpha}\right)$$
+
+$$\text{FD}_{\alpha\beta} = N_\alpha N_\beta \sin\frac{2\pi m_\alpha}{N_\alpha} \sin\frac{2\pi m_\beta}{N_\beta} \quad (\alpha \neq \beta)$$
+
+FD Laplacian 的等效 gg 为：
+
+$$\text{gg}_{\text{FD}} = \frac{1}{(2\pi)^2} \sum_{\alpha,\beta} \text{GGT}[\alpha][\beta] \cdot \text{FD}_{\alpha\beta}$$
+
+最终：
+
+$$\nabla^2_{\text{FD}} f = -\text{gg}_{\text{FD}} \cdot \text{tpiba2} \cdot f$$
+
+**与谱 kernel 的关系**：当 $m_\alpha / N_\alpha \to 0$ 时，
+
+$$2(1-\cos\theta) \to \theta^2, \quad \sin\theta \to \theta$$
+
+因此 $\text{FD}_{\alpha\alpha} \to (2\pi m_\alpha)^2$，$\text{FD}_{\alpha\beta} \to (2\pi)^2 m_\alpha m_\beta$，$\text{gg}_{\text{FD}} \to \text{gg}_{\text{spectral}}$。
+
+**有界性**：gg_FD 的最大值有限（由 FFT 网格大小决定），而 gg_spectral 随 |m| 无界增长。
+
+### 2.6 正交晶胞简化
+
+对正交晶胞，GGT 为对角矩阵（GGT[αβ] = 0 for α ≠ β），交叉项消失：
+
+$$\text{gg}_{\text{FD}} = \frac{1}{(2\pi)^2} \sum_\alpha \text{GGT}[\alpha][\alpha] \cdot 2N_\alpha^2\left(1 - \cos\frac{2\pi m_\alpha}{N_\alpha}\right)$$
+
+每个方向独立贡献。
+
+---
+
+## 3. 实现
+
+### 3.1 vlapl 势
+
+在 `v_xc_meta` 函数中，对每个自旋通道 is：
+
+1. **构造 vlapl·sgn**（实空间）：
+
+$$\text{vlapl\_sgn}[ir] = \text{vlapl}[ir \cdot \text{nspin} + \text{is}] \cdot \text{sgn}[ir \cdot \text{nspin} + \text{is}]$$
+
+2. **FFT 到 G 空间**：
+
+$$\widetilde{\text{vlapl\_sgn}} = \text{FFT}(\text{vlapl\_sgn})$$
+
+3. **应用 FD kernel**：对每个 G 向量 ig，Miller 指数 m = gdirect[ig]，
+
+$$\widetilde{\nabla^2_{\text{FD}}(\text{vlapl\_sgn})}[ig] = -\text{gg}_{\text{FD}}[ig] \cdot \text{tpiba2} \cdot \widetilde{\text{vlapl\_sgn}}[ig]$$
+
+其中 gg_FD[ig] 由 2.5 节公式计算。
+
+4. **IFFT 回实空间**：
+
+$$\nabla^2_{\text{FD}}(\text{vlapl\_sgn}) = \text{IFFT}(\widetilde{\nabla^2_{\text{FD}}(\text{vlapl\_sgn})})$$
+
+5. **加入 V_xc**：
+
+$$V_{xc}(\text{is}, ir) \mathrel{+}= e^2 \cdot \nabla^2_{\text{FD}}(\text{vlapl\_sgn})[ir]$$
+
+对于 HSE 等 hybrid 泛函 (func_type == 5)，乘以 (1 − hybrid_alpha)。
+
+6. **加入 vtxc**：
+
+$$\text{vtxc} \mathrel{+}= e^2 \sum_{ir} \nabla^2_{\text{FD}}(\text{vlapl\_sgn})[ir] \cdot \rho[\text{is}][ir] \cdot \text{sgn}[ir \cdot \text{nspin} + \text{is}]$$
+
+### 3.2 vlapl 应力
+
+vlapl 对应力的贡献在 `gradcorr` 函数中计算，与谱 Laplacian 方案相同：
+
+$$\sigma_{\alpha\beta}^{\nabla^2\rho} = \frac{2}{\Omega} \sum_{\mathbf{r}} \text{vlapl}(\mathbf{r}) \cdot H_{\alpha\beta}(\mathbf{r}) \cdot e^2$$
+
+其中 $H_{\alpha\beta} = \partial^2\rho / \partial r_\alpha \partial r_\beta$ 是密度的 Hessian 矩阵（谱方法计算）。LibXC 返回的 vlapl = ∂(ρε)/∂(∇²ρ) 已包含 ρ 因子。
+
+### 3.3 功能名注册
+
+`dft_functional scanl` 映射为 LibXC 泛函 `XC_MGGA_X_SCANL` (ID=700) + `XC_MGGA_C_SCANL` (ID=702)。
+
+---
+
+## 4. FD kernel 数值验证
+
+### 4.1 验证设置
+
+FCC 原胞 (Si₂)，a = 10.2 Bohr。
+
+晶格矢量（lat0 单位）：
+
+$$\mathbf{a}_1 = (0, \tfrac{1}{2}, \tfrac{1}{2}), \quad \mathbf{a}_2 = (\tfrac{1}{2}, 0, \tfrac{1}{2}), \quad \mathbf{a}_3 = (\tfrac{1}{2}, \tfrac{1}{2}, 0)$$
+
+度规矩阵：
+
+$$\text{GGT} = \begin{pmatrix} 3 & -1 & -1 \\ -1 & 3 & -1 \\ -1 & -1 & 3 \end{pmatrix}$$
+
+### 4.2 gg_FD vs gg_spectral（网格 36×36×36, ecutwfc=60 Ry）
+
+| Miller (m₁,m₂,m₃) | gg_spectral | gg_FD | gg_FD / gg_spectral |
+|---------------------|-------------|-------|---------------------|
+| (0,0,0) | 0 | 0 | — |
+| (1,0,0) | 3.000 | 2.992 | 0.9975 |
+| (1,1,0) | 4.000 | 4.005 | 1.0013 |
+| (2,0,0) | 12.00 | 11.88 | 0.9900 |
+| (3,0,0) | 27.00 | 26.46 | 0.9800 |
+| (5,0,0) | 75.00 | 70.36 | 0.9381 |
+| (8,0,0) | 192.0 | 163.4 | 0.8511 |
+| (10,0,0) | 300.0 | 231.2 | 0.7706 |
+| (15,0,0) | 675.0 | 387.4 | 0.5740 |
+| (18,0,0) | 972.0 | 393.9 | 0.4053 |
+
+**解读**：
+
+- |m|/N ≤ 0.1 时：FD 与谱偏差 < 1%，物理上等价
+- |m|/N ≈ 0.2 时：偏差 ~6%，开始有差异
+- |m|/N ≈ 0.5（Nyquist）：FD = 0.405 × 谱，大幅衰减
+
+### 4.3 gg_FD vs gg_spectral（网格 45×45×45, ecutwfc=80 Ry）
+
+| Miller (m₁,m₂,m₃) | gg_spectral | gg_FD | gg_FD / gg_spectral |
+|---------------------|-------------|-------|---------------------|
+| (1,0,0) | 3.000 | 2.994 | 0.9980 |
+| (3,0,0) | 27.00 | 26.66 | 0.9875 |
+| (5,0,0) | 75.00 | 72.55 | 0.9673 |
+| (8,0,0) | 192.0 | 173.7 | 0.9048 |
+| (9,0,0) | 243.0 | 211.7 | 0.8712 |
+
+ecutwfc=80 对应的最大 |m| ≈ 8.4，此时 FD/谱 ≈ 0.89。即 PW 截断处 FD kernel 衰减约 11%。
+
+### 4.4 有界性验证
+
+Nyquist 频率 (m_α = N_α/2) 下 FD 因子的最大值：
+
+$$\text{FD}_{\alpha\alpha}^{\max} = 2N_\alpha^2 \cdot 2 = 4N_\alpha^2$$
+
+$$\text{gg}_{\text{FD}}^{\max} = \frac{1}{(2\pi)^2} \sum_\alpha \text{GGT}[\alpha][\alpha] \cdot 4N_\alpha^2 = \frac{4}{(2\pi)^2} \cdot 3 \cdot 3 N^2 = \frac{36 N^2}{4\pi^2}$$
+
+对 N=36：gg_FD_max = 1181.8，对应 |G|_FD_max ≈ 21.2 Bohr⁻¹（有界）。
+对 N=45：gg_FD_max = 1846.5。
+
+而谱 Laplacian 在同一点：gg_spectral = 9 · (N/2)² → 随 N² 无界增长。
+
+---
+
+## 5. 有限差分应力验证
+
+### 5.1 验证方法
+
+对 Si₂ FCC 原胞（Γ 点）在 a₀ ± δ 下计算总能量 E，由差分得到应力：
+
+$$\sigma_{\text{FD}} = -\frac{E(a_0+\delta) - E(a_0-\delta)}{V(a_0+\delta) - V(a_0-\delta)}$$
+
+FCC 原胞体积 $V = |\det(\text{latvec})| \cdot a^3 / 4 = a^3 / 4$。
+
+与 ABACUS 解析应力 σ_AB 比较。单位转换：
+
+$$\sigma(\text{kbar}) = \sigma(\text{eV/Bohr}^3) \times \frac{147105}{13.6057}$$
+
+其中 147105 kbar/(Ry/Bohr³)，13.6057 eV/Ry。
+
+### 5.2 测试参数
+
+| 参数 | 值 |
+|------|-----|
+| 体系 | Si₂ FCC 原胞 |
+| a₀ | 10.2 Bohr |
+| δ | 0.001 Bohr |
+| 赝势 | Si_ONCV_PBE-1.0 (NCPP) |
+| smearing | gauss, σ = 0.002 |
+| mixing | Pulay, beta = 0.1–0.2, gg0 = 1.0 |
+| mixing_tau | 1 |
+| scf_thr | 1e-8 |
+| pw_seed | 1 |
+
+### 5.3 SCAN 基线（无 ∇²ρ 依赖）
+
+SCAN 的 vlapl = 0，FD Laplacian 代码路径不激活，用于验证 FD 方法和应力代码。
+
+| ecutwfc (Ry) | FFT grid | σ_FD (kbar) | σ_AB (kbar) | FD/AB | 误差 |
+|-------------|----------|-------------|-------------|-------|------|
+| 60 | 36³ | 434.29 | 434.30 | 1.0000 | 0.00% |
+| 80 | 45³ | 434.47 | 434.48 | 1.0000 | 0.00% |
+| 100 | 45³ | 434.54 | 434.55 | 1.0000 | 0.00% |
+
+**结论**：SCAN 在所有 ecutwfc 下 FD/AB = 1.0000，验证了代码和方法的正确性。
+
+### 5.4 SCANL（FD Laplacian 方案）
+
+| ecutwfc (Ry) | FFT grid | σ_FD (kbar) | σ_AB (kbar) | FD/AB | 误差 |
+|-------------|----------|-------------|-------------|-------|------|
+| 60 | 36³ | 431.89 | 435.60 | 0.9915 | -0.85% |
+| 80 | 45³ | 432.37 | 434.58 | 0.9949 | -0.51% |
+
+ecutwfc=100 时 SCANL 的 SCF 收敛困难（能量振荡），FD 结果不可靠，未列入。
+
+### 5.5 完整能量数据
+
+ecutwfc = 80 Ry：
+
+| | E(a₀+δ) (eV) | E(a₀−δ) (eV) | ΔE (eV) |
+|------|-------------|-------------|---------|
+| SCAN | −197.22465 | −197.21838 | −0.006271 |
+| SCANL | −198.70596 | −198.69972 | −0.006241 |
+
+ecutwfc = 60 Ry：
+
+| | E(a₀+δ) (eV) | E(a₀−δ) (eV) | ΔE (eV) |
+|------|-------------|-------------|---------|
+| SCAN | −197.22344 | −197.21717 | −0.006269 |
+| SCANL | −198.68669 | −198.68046 | −0.006234 |
+
+SCAN 的 ΔE 跨 ecutwfc 稳定（−0.00627 eV），SCANL 的 ΔE 在 ecut=60/80 时接近（−0.00623/−0.00624 eV），与 SCAN 的差异约 0.5%，与 FD 应力误差一致。
+
+### 5.6 与其他方案对比
+
+| 方案 | V 含 vlapl | vtxc 含 vlapl | stress 含 vlapl | SCF | FD/AB (ecut=80) | 误差 |
+|------|-----------|--------------|----------------|-----|-------|------|
+| 严格谱 Laplacian | ✓ (谱) | ✓ (谱) | ✓ | **发散** | — | — |
+| sigma-as-lapl | ✗ (假输入) | ✗ (假输入) | ✗ | 收敛 | 1.046 | +4.6% |
+| V/vtxc 均省略 vlapl | ✗ | ✗ | ✓ | 收敛 | 0.906 | −9.5% |
+| vtxc 含 vlapl | ✗ | ✓ | ✓ | 收敛 | 0.922 | −7.8% |
+| 高斯衰减 | ✓ (衰减) | ✓ (衰减) | ✓ | 收敛 | ~0.9 | ~10% |
+| 绝热 ramp | ✓ (渐增) | ✓ | ✓ | **发散** | — | — |
+| **FD Laplacian** | **✓ (FD)** | **✓ (FD)** | **✓** | **收敛** | **0.995** | **−0.5%** |
+
+---
+
+## 6. 误差分析
+
+### 6.1 残余 0.5% 应力误差来源
+
+1. **FD kernel 与谱 kernel 在高 G 的偏差**：在 PW 截断处（|m|/N ≈ 0.37 for ecut=80），FD kernel 衰减到谱值的 ~89%。vlapl 势在高 G 分量上偏弱，密度与严格 SCANL 密度有微小差异。
+
+2. **FD 应力本身的高阶误差**：δ = 0.001 Bohr 的中心差分为二阶精度，误差 O(δ²/a⁴) ≈ 10⁻⁸，对 0.5% 误差贡献可忽略。
+
+3. **vlapl stress 项中 Hessian 仍为谱 Hessian**：gradcorr 中的 vlapl 应力项使用谱方法计算 Hessian，而势用 FD Laplacian。此不一致性贡献未知但量级较小。
+
+### 6.2 ecutwfc=100 的 SCF 困难
+
+ecutwfc=100 时 SCANL 的 ΔE = −0.00568 eV，显著偏离 SCAN 的 −0.00628 eV（偏差 ~10%），而 ecut=60/80 时偏差仅 ~0.5%。原因可能是：
+
+- 45³ FFT 网格对 ecut=100 的 35749 个平面波不够大，某些高 G 分量的 FD kernel 在网格边界处行为异常
+- SCF 收敛到与 ecut=60/80 不同的密度
+
+**建议**：SCANL 使用 ecutwfc ≤ 80 Ry。若需更高截断，应确保 FFT 网格足够大。
+
+### 6.3 与严格实现的关系
+
+严格 SCANL 实现要求谱 Laplacian ∇² = −|G|²，但在任何实际 ecutwfc 下导致 SCF 发散。FD Laplacian 是在 FFT 网格上对离散函数求 Laplacian 的自洽方法——它是网格函数的"正确" Laplacian，但与连续 Laplacian 有差异。此差异随网格加密而减小（FD → spectral 当 N → ∞），但在实际网格上不可消除。
+
+---
+
+## 7. 使用建议
+
+1. **dft_functional scanl**：已注册，映射到 XC_MGGA_X_SCANL + XC_MGGA_C_SCANL
+2. **ecutwfc ≤ 80 Ry**：ecutwfc 过高时 SCF 收敛困难
+3. **mixing_tau = 1**：meta-GGA 必需
+4. **Pulay mixing, beta = 0.1–0.2**：对 SCANL 收敛效果较好
+5. **优先考虑 r²SCAN 或 SCAN**：无 ∇²ρ 依赖，数值更稳定
diff --git a/source/module_hamilt_general/module_xc/test/CMakeLists.txt b/source/module_hamilt_general/module_xc/test/CMakeLists.txt
index 0dda934ac6b..5cb216d5e54 100644
--- a/source/module_hamilt_general/module_xc/test/CMakeLists.txt
+++ b/source/module_hamilt_general/module_xc/test/CMakeLists.txt
@@ -11,6 +11,7 @@ AddTest(
   TARGET XCTest_HSE
   LIBS parameter MPI::MPI_CXX Libxc::xc  # required by global.h; for details, `remove_definitions(-D__MPI)`.
   SOURCES test_xc1.cpp ../xc_functional.cpp ../xc_functional_libxc.cpp
+    ../xc_functional_libxc_tools.cpp
 )
 
 
@@ -36,6 +37,7 @@ AddTest(
     ../xc_functional_libxc_wrapper_xc.cpp
     ../xc_functional_libxc_wrapper_gcxc.cpp
     ../xc_functional_libxc_wrapper_tauxc.cpp
+    ../xc_functional_libxc_tools.cpp
     ../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp ../xc_funct_exch_gga.cpp
     ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp
     ../../../module_base/matrix.cpp
@@ -57,6 +59,7 @@ AddTest(
     ../xc_functional_libxc_wrapper_xc.cpp
     ../xc_functional_libxc_wrapper_gcxc.cpp
     ../xc_functional_libxc_wrapper_tauxc.cpp
+    ../xc_functional_libxc_tools.cpp
     ../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp 
     ../xc_funct_exch_gga.cpp ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp 
 )
diff --git a/source/module_hamilt_general/module_xc/test/test_xc4.cpp b/source/module_hamilt_general/module_xc/test/test_xc4.cpp
index 41fd60f4254..113bae779f3 100644
--- a/source/module_hamilt_general/module_xc/test/test_xc4.cpp
+++ b/source/module_hamilt_general/module_xc/test/test_xc4.cpp
@@ -43,8 +43,9 @@ class XCTest_SCAN : public XCTest
 
             for(int i=0;i<5;i++)
             {
-                double e,v,v1,v2,v3;
-                XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],tau[i],e,v1,v2,v3);
+                double e,v,v1,v2,v3,v4;
+                // SCAN does not depend on laplacian, pass 0.0
+                XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],0.0,tau[i],e,v1,v2,v3,v4);
                 e_.push_back(e);
                 v1_.push_back(v1);
                 v2_.push_back(v2);
diff --git a/source/module_hamilt_general/module_xc/xc_functional.cpp b/source/module_hamilt_general/module_xc/xc_functional.cpp
index 2fc74ffbb78..e8758916d4f 100644
--- a/source/module_hamilt_general/module_xc/xc_functional.cpp
+++ b/source/module_hamilt_general/module_xc/xc_functional.cpp
@@ -199,6 +199,13 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
         func_type = 5;
         use_libxc = true;
     }
+    else if ( xc_func == "SCANL")
+    {
+        func_id.push_back(XC_MGGA_X_SCANL);
+        func_id.push_back(XC_MGGA_C_SCANL);
+        func_type = 3;
+        use_libxc = true;
+    }
     else if( xc_func == "LC_PBE")
     {
         func_id.push_back(XC_HYB_GGA_XC_LC_PBEOP);
@@ -350,7 +357,7 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
 
 #ifndef USE_LIBXC
     if(xc_func == "SCAN" || xc_func == "HSE" || xc_func == "SCAN0" 
-        || xc_func == "MULLER" || xc_func == "POWER" || xc_func == "WP22" || xc_func == "CWP22" ||
+        || xc_func == "SCANL" || xc_func == "MULLER" || xc_func == "POWER" || xc_func == "WP22" || xc_func == "CWP22" ||
         xc_func == "LC_PBE" || xc_func == "LC_WPBE" || xc_func == "LRC_WPBE" ||
         xc_func == "LRC_PBEH" || xc_func == "CAM_PBEH")
     {
diff --git a/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp b/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp
index 42501465108..b1a8f9501a9 100644
--- a/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp
+++ b/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp
@@ -67,6 +67,10 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 	ModuleBase::Vector3<double>* gdr2 = nullptr;
 	ModuleBase::Vector3<double>* h1 = nullptr;
 	ModuleBase::Vector3<double>* h2 = nullptr;
+	double* lapl1 = nullptr;
+	double* lapl2 = nullptr;
+	double* hess1 = nullptr; // nspin=1 Hessian (6 components)
+	double* hess2 = nullptr; // nspin=1 Hessian for spin-down channel
 	double* neg = nullptr;
 	double** vsave = nullptr;
 	double** vgg = nullptr;
@@ -95,6 +99,46 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 	
 	XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba);
 
+	XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba);
+
+	// compute laplacian of total density (rho[0] + rho_core) for nspin=1,2
+	// and as the "up" channel equivalent for nspin=4 noncollinear case
+	{
+		std::vector<double> rho_total(rhopw->nrxx);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+		for(int ir=0; ir<rhopw->nrxx; ++ir)
+			rho_total[ir] = rhotmp1[ir];
+		std::vector<double> lapl_all = XC_Functional_Libxc::cal_lapl(
+			1, rhopw->nrxx, rho_total, ucell->tpiba2, chr);
+		lapl1 = new double[rhopw->nrxx];
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+		for(int ir=0; ir<rhopw->nrxx; ++ir)
+			lapl1[ir] = lapl_all[ir];
+	}
+
+	// compute Hessian of total density for meta-GGA stress
+	if(is_stress && (func_type == 3 || func_type == 5))
+	{
+		std::vector<double> rho_total(rhopw->nrxx);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+		for(int ir=0; ir<rhopw->nrxx; ++ir)
+			rho_total[ir] = rhotmp1[ir];
+		std::vector<double> hess_all = XC_Functional_Libxc::cal_rho_hessian(
+			1, rhopw->nrxx, rho_total, chr);
+		hess1 = new double[rhopw->nrxx * 6];
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+		for(int i=0; i<rhopw->nrxx * 6; ++i)
+			hess1[i] = hess_all[i];
+	}
+
 	// for spin polarized case;
 	// calculate the gradient of (rho_core+rho) in reciprocal space.
 	if(PARAM.inp.nspin==2)
@@ -120,6 +164,43 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 		if(!is_stress) { h2 = new ModuleBase::Vector3<double>[rhopw->nrxx]; }
 		
 		XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba);
+
+		// compute laplacian of spin-down density (rho[1] + rho_core)
+		{
+			std::vector<double> rho_dw(rhopw->nrxx);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int ir=0; ir<rhopw->nrxx; ++ir)
+				rho_dw[ir] = rhotmp2[ir];
+			std::vector<double> lapl_all = XC_Functional_Libxc::cal_lapl(
+				1, rhopw->nrxx, rho_dw, ucell->tpiba2, chr);
+			lapl2 = new double[rhopw->nrxx];
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int ir=0; ir<rhopw->nrxx; ++ir)
+				lapl2[ir] = lapl_all[ir];
+		}
+
+		// compute Hessian of spin-down density for meta-GGA stress
+		if(is_stress && (func_type == 3 || func_type == 5))
+		{
+			std::vector<double> rho_dw(rhopw->nrxx);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int ir=0; ir<rhopw->nrxx; ++ir)
+				rho_dw[ir] = rhotmp2[ir];
+			std::vector<double> hess_all = XC_Functional_Libxc::cal_rho_hessian(
+				1, rhopw->nrxx, rho_dw, chr);
+			hess2 = new double[rhopw->nrxx * 6];
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int i=0; i<rhopw->nrxx * 6; ++i)
+				hess2[i] = hess_all[i];
+		}
 	}
 
 	if(PARAM.inp.nspin == 4&&(PARAM.globalv.domag||PARAM.globalv.domag_z))
@@ -189,6 +270,40 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 		XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba);
 		XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba);
 
+		// re-compute laplacians after noncollinear rho update
+		{
+			std::vector<double> rho_up(rhopw->nrxx);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int ir=0; ir<rhopw->nrxx; ++ir)
+				rho_up[ir] = rhotmp1[ir];
+			std::vector<double> lapl_up = XC_Functional_Libxc::cal_lapl(
+				1, rhopw->nrxx, rho_up, ucell->tpiba2, chr);
+			lapl1 = new double[rhopw->nrxx];
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int ir=0; ir<rhopw->nrxx; ++ir)
+				lapl1[ir] = lapl_up[ir];
+		}
+		{
+			std::vector<double> rho_dw(rhopw->nrxx);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int ir=0; ir<rhopw->nrxx; ++ir)
+				rho_dw[ir] = rhotmp2[ir];
+			std::vector<double> lapl_dw = XC_Functional_Libxc::cal_lapl(
+				1, rhopw->nrxx, rho_dw, ucell->tpiba2, chr);
+			lapl2 = new double[rhopw->nrxx];
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+			for(int ir=0; ir<rhopw->nrxx; ++ir)
+				lapl2[ir] = lapl_dw[ir];
+		}
+
 	}
 	
 	const double epsr = 1.0e-6;
@@ -249,9 +364,20 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 #ifdef USE_LIBXC
 					if(func_type == 3 || func_type == 5) //the gradcorr part to stress of mGGA
 					{
-						double v3xc;
+						double v3xc, vlaplc;
 						double atau = chr->kin_r[0][ir]/2.0;
-						XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, atau, sxc, v1xc, v2xc, v3xc);
+						XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, lapl1[ir], atau, sxc, v1xc, v2xc, v3xc, vlaplc);
+						// vlapl stress: σ_αβ = (2/Ω) Σ vlapl·Hessian_αβ·e2
+						// vlapl = ∂(ρε)/∂(∇²ρ) already includes ρ factor
+						for(int l = 0; l < 3; l++)
+						{
+							for(int m = 0; m < l+1; m++)
+							{
+								int ind = l*3 + m;
+								int ic = (l==0&&m==0) ? 0 : (l==1&&m==1) ? 1 : (l==2&&m==2) ? 2 : (l==0&&m==1) ? 3 : (l==1&&m==2) ? 4 : 5;
+								local_stress_gga[ind] += 2.0 * vlaplc * hess1[ic * rhopw->nrxx + ir] * ModuleBase::e2;
+							}
+						}
 					}
 					else
 					{
@@ -308,13 +434,28 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 				double sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud;
 				if(func_type == 3 || func_type == 5) //the gradcorr part to stress of mGGA
 				{
-					double v3xcup, v3xcdw;
+					double v3xcup, v3xcdw, vlaplup, vlapldw;
 					double atau1 = chr->kin_r[0][ir]/2.0;
 					double atau2 = chr->kin_r[1][ir]/2.0;
 					XC_Functional_Libxc::tau_xc_spin(
 						func_id,
-						rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], 
-						atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw);
+						rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir],
+						lapl1[ir], lapl2[ir],
+						atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, vlaplup, vlapldw);
+					// vlapl stress: σ_αβ = (2/Ω) Σ (vlapl_up·Hess_up + vlapl_dw·Hess_dw)·e2
+					// vlapl = ∂(ρε)/∂(∇²ρ) already includes ρ factor
+					if(is_stress)
+					{
+						for(int l = 0; l < 3; l++)
+						{
+							for(int m = 0; m < l+1; m++)
+							{
+								int ind = l*3 + m;
+								int ic = (l==0&&m==0) ? 0 : (l==1&&m==1) ? 1 : (l==2&&m==2) ? 2 : (l==0&&m==1) ? 3 : (l==1&&m==2) ? 4 : 5;
+								local_stress_gga[ind] += 2.0 * (vlaplup * hess1[ic * rhopw->nrxx + ir] + vlapldw * hess2[ic * rhopw->nrxx + ir]) * ModuleBase::e2;
+							}
+						}
+					}
 				}
 				else
 				{
@@ -582,6 +723,8 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 	delete[] rhotmp1;
 	delete[] rhogsum1;
 	delete[] gdr1;
+	delete[] lapl1;
+	delete[] hess1;
 	if(!is_stress) { delete[] h1;
 }
 
@@ -590,6 +733,8 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 		delete[] rhotmp2;
 		delete[] rhogsum2;
 		delete[] gdr2;
+		delete[] lapl2;
+		delete[] hess2;
 		if(!is_stress) { delete[] h2;
 }
 	}
@@ -609,6 +754,8 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
 		delete[] rhotmp2;
 		delete[] rhogsum2;
 		delete[] gdr2;
+		delete[] lapl2;
+		delete[] hess2;
 	}
 
 	return;
diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc.h b/source/module_hamilt_general/module_xc/xc_functional_libxc.h
index 12cf1fdd8ed..acadcd866f9 100644
--- a/source/module_hamilt_general/module_xc/xc_functional_libxc.h
+++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.h
@@ -82,6 +82,22 @@ namespace XC_Functional_Libxc
     extern std::vector<double> convert_sigma(
         const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr);
 
+    // calculating laplacian of density: ∇²ρ(r)
+    extern std::vector<double> cal_lapl(
+        const int nspin,
+        const std::size_t nrxx,
+        const std::vector<double> &rho,
+        const double tpiba2,
+        const Charge* const chr);
+
+    // calculating Hessian of density: ∂²ρ/∂r_α∂r_β
+    // returns nspin * nrxx * 6 array (xx, yy, zz, xy, yz, zx)
+    extern std::vector<double> cal_rho_hessian(
+        const int nspin,
+        const std::size_t nrxx,
+        const std::vector<double> &rho,
+        const Charge* const chr);
+
     // sgn for threshold mask
     extern std::vector<double> cal_sgn(
         const double rho_threshold,
@@ -166,16 +182,17 @@ namespace XC_Functional_Libxc
     // wrapper for the mGGA functionals
     extern void tau_xc(
         const std::vector<int> &func_id,
-        const double &rho, const double &grho, const double &atau, double &sxc,
-        double &v1xc, double &v2xc, double &v3xc);
+        const double &rho, const double &grho, const double &lapl_rho, const double &atau, double &sxc,
+        double &v1xc, double &v2xc, double &v3xc, double &vlaplc);
 
     extern void tau_xc_spin(
         const std::vector<int> &func_id,
         double rhoup, double rhodw,
         ModuleBase::Vector3<double> gdr1, ModuleBase::Vector3<double> gdr2,
+        double laplup, double lapldw,
         double tauup, double taudw,
         double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud,
-        double &v3xcup, double &v3xcdw);
+        double &v3xcup, double &v3xcdw, double &vlaplup, double &vlapldw);
 } // namespace XC_Functional_Libxc
 
 #endif // USE_LIBXC
diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp
index 567d04e62a8..f4297635cfd 100644
--- a/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp
+++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp
@@ -82,6 +82,98 @@ XC_Functional_Libxc::cal_gdr(
 	return gdr;
 }
 
+// calculating laplacian of density: ∇²ρ(r)
+std::vector<double> XC_Functional_Libxc::cal_lapl(
+	const int nspin,
+	const std::size_t nrxx,
+	const std::vector<double> &rho,
+	const double tpiba2,
+	const Charge* const chr)
+{
+	std::vector<double> lapl(nspin * nrxx, 0.0);
+	for( int is=0; is!=nspin; ++is )
+	{
+		std::vector<double> rhor(nrxx);
+		#ifdef _OPENMP
+		#pragma omp parallel for schedule(static, 1024)
+		#endif
+		for(std::size_t ir=0; ir<nrxx; ++ir)
+			rhor[ir] = rho[ir*nspin+is];
+
+		std::vector<std::complex<double>> rhog(chr->rhopw->npw);
+		chr->rhopw->real2recip(rhor.data(), rhog.data());
+
+		// ∇²ρ(r) = -∑_G |G|² ρ(G) e^{iGr}
+		std::vector<std::complex<double>> rhog_lapl(chr->rhopw->npw);
+		#ifdef _OPENMP
+		#pragma omp parallel for schedule(static, 1024)
+		#endif
+		for(int ig=0; ig<chr->rhopw->npw; ++ig)
+			rhog_lapl[ig] = -rhog[ig] * chr->rhopw->gg[ig] * tpiba2;
+
+		std::vector<std::complex<double>> aux(chr->rhopw->nmaxgr);
+		chr->rhopw->recip2real(rhog_lapl.data(), aux.data());
+		#ifdef _OPENMP
+		#pragma omp parallel for schedule(static, 1024)
+		#endif
+		for(std::size_t ir=0; ir<nrxx; ++ir)
+			lapl[is * nrxx + ir] = aux[ir].real();
+	} // end for(is)
+	return lapl;
+}
+
+// calculating Hessian of density: ∂²ρ/∂r_α∂r_β
+// returns nspin * nrxx * 6 array (xx, yy, zz, xy, yz, zx)
+std::vector<double> XC_Functional_Libxc::cal_rho_hessian(
+	const int nspin,
+	const std::size_t nrxx,
+	const std::vector<double> &rho,
+	const Charge* const chr)
+{
+	std::vector<double> hess(nspin * nrxx * 6, 0.0);
+	// ipol2xy maps 6 components to (x,y) pairs: 0:xx, 1:yy, 2:zz, 3:xy, 4:yz, 5:zx
+	const int ipol2xy[6][2] = {{0,0}, {1,1}, {2,2}, {0,1}, {1,2}, {2,0}};
+	
+	for( int is=0; is!=nspin; ++is )
+	{
+		std::vector<double> rhor(nrxx);
+		#ifdef _OPENMP
+		#pragma omp parallel for schedule(static, 1024)
+		#endif
+		for(std::size_t ir=0; ir<nrxx; ++ir)
+			rhor[ir] = rho[ir*nspin+is];
+
+		std::vector<std::complex<double>> rhog(chr->rhopw->npw);
+		chr->rhopw->real2recip(rhor.data(), rhog.data());
+
+		// compute all 6 Hessian components in one G-space pass
+		for(int ic=0; ic<6; ++ic)
+		{
+			const int ax = ipol2xy[ic][0];
+			const int ay = ipol2xy[ic][1];
+			std::vector<std::complex<double>> rhog_hess(chr->rhopw->npw);
+			#ifdef _OPENMP
+			#pragma omp parallel for schedule(static, 1024)
+			#endif
+			for(int ig=0; ig<chr->rhopw->npw; ++ig)
+			{
+				const double gx = chr->rhopw->gcar[ig][ax];
+				const double gy = chr->rhopw->gcar[ig][ay];
+				rhog_hess[ig] = -rhog[ig] * gx * gy;
+			}
+
+			std::vector<std::complex<double>> aux(chr->rhopw->nmaxgr);
+			chr->rhopw->recip2real(rhog_hess.data(), aux.data());
+			#ifdef _OPENMP
+			#pragma omp parallel for schedule(static, 1024)
+			#endif
+			for(std::size_t ir=0; ir<nrxx; ++ir)
+				hess[is * nrxx * 6 + ic * nrxx + ir] = aux[ir].real();
+		}
+	} // end for(is)
+	return hess;
+}
+
 // converting grho (abacus=>libxc)
 std::vector<double> XC_Functional_Libxc::convert_sigma(
 	const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr)
diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp
index 2456abe5ba9..6c6558615bd 100644
--- a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp
+++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp
@@ -199,6 +199,10 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
         = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr);
     const std::vector<double> sigma = XC_Functional_Libxc::convert_sigma(gdr);
 
+    // compute laplacian of density: ∇²ρ(r) for each spin channel
+    const double tpiba2 = tpiba * tpiba;
+    const std::vector<double> lapl = XC_Functional_Libxc::cal_lapl(nspin, nrxx, rho, tpiba2, chr);
+
     //converting kin_r
     std::vector<double> kin_r;
     kin_r.resize(nrxx*nspin);
@@ -264,7 +268,7 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
     for ( xc_func_type &func : funcs )
     {
         assert(func.info->family == XC_FAMILY_MGGA);
-        xc_mgga_exc_vxc(&func, nrxx, rho.data(), sigma.data(), sigma.data(),
+        xc_mgga_exc_vxc(&func, nrxx, rho.data(), sigma.data(), lapl.data(),
             kin_r.data(), exc.data(), vrho.data(), vsigma.data(), vlapl.data(), vtau.data());
 
         //process etxc
@@ -387,6 +391,99 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
                 vofk(is,ir) += vtau[ir*nspin+is]  * sgn[ir*nspin+is];
             }
         }
+
+        //process vlapl: ∇²(vlapl) contribution to potential
+        // Use finite-difference Laplacian in G-space to avoid unbounded
+        // |G|² amplification from the spectral Laplacian. The FD kernel
+        // replaces m_α² with n_α²·2(1-cos(2πm_α/n_α)) and m_α·m_β
+        // with n_α·n_β·sin(2πm_α/n_α)·sin(2πm_β/n_β), which matches
+        // the spectral kernel at low G but is bounded at high G.
+        // Cross terms from the metric tensor (GGT) are included for
+        // non-orthogonal cells.
+        for( int is=0; is<nspin; ++is )
+        {
+            std::vector<double> vlapl_sgn(nrxx);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+            for( int ir=0; ir<nrxx; ++ir )
+            {
+                vlapl_sgn[ir] = vlapl[ir*nspin+is] * sgn[ir*nspin+is];
+            }
+
+            std::vector<std::complex<double>> vlapl_g(chr->rhopw->npw);
+            chr->rhopw->real2recip(vlapl_sgn.data(), vlapl_g.data());
+
+            const int ngrid[3] = {chr->rhopw->nx, chr->rhopw->ny, chr->rhopw->nz};
+            const ModuleBase::Matrix3& GGT_mat = chr->rhopw->GGT;
+            const double two_pi = ModuleBase::TWO_PI;
+            const double ggt[3][3] = {
+                {GGT_mat.e11, GGT_mat.e12, GGT_mat.e13},
+                {GGT_mat.e21, GGT_mat.e22, GGT_mat.e23},
+                {GGT_mat.e31, GGT_mat.e32, GGT_mat.e33}
+            };
+
+            const double four_pi_sq = ModuleBase::FOUR_PI * ModuleBase::PI;
+            std::vector<std::complex<double>> vlapl_g_lapl(chr->rhopw->npw);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 256)
+#endif
+            for(int ig=0; ig<chr->rhopw->npw; ++ig)
+            {
+                double gg_fd_raw = 0.0;
+                const double m[3] = {chr->rhopw->gdirect[ig].x,
+                                     chr->rhopw->gdirect[ig].y,
+                                     chr->rhopw->gdirect[ig].z};
+                for(int alpha=0; alpha<3; ++alpha)
+                {
+                    const double phi_a = two_pi * m[alpha] / ngrid[alpha];
+                    for(int beta=0; beta<3; ++beta)
+                    {
+                        const double phi_b = two_pi * m[beta] / ngrid[beta];
+                        if(alpha == beta)
+                        {
+                            gg_fd_raw += ggt[alpha][alpha] * ngrid[alpha] * ngrid[alpha] * 2.0 * (1.0 - cos(phi_a));
+                        }
+                        else
+                        {
+                            gg_fd_raw += ggt[alpha][beta] * ngrid[alpha] * ngrid[beta] * sin(phi_a) * sin(phi_b);
+                        }
+                    }
+                }
+                const double gg_fd = gg_fd_raw / four_pi_sq;
+                vlapl_g_lapl[ig] = -vlapl_g[ig] * gg_fd * tpiba2;
+            }
+
+            std::vector<std::complex<double>> aux(chr->rhopw->nmaxgr);
+            chr->rhopw->recip2real(vlapl_g_lapl.data(), aux.data());
+
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 1024)
+#endif
+            for( int ir=0; ir<nrxx; ++ir )
+            {
+#ifdef __EXX
+                double vlapl_contrib = aux[ir].real();
+                if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5)
+                {
+                    vlapl_contrib *= (1.0 - XC_Functional::get_hybrid_alpha());
+                }
+                v(is,ir) += ModuleBase::e2 * vlapl_contrib;
+#else
+                v(is,ir) += ModuleBase::e2 * aux[ir].real();
+#endif
+            }
+
+            double vtxc_vlapl = 0.0;
+#ifdef _OPENMP
+#pragma omp parallel for reduction(+:vtxc_vlapl) schedule(static, 1024)
+#endif
+            for( int ir=0; ir<nrxx; ++ir )
+            {
+                vtxc_vlapl += ModuleBase::e2 * aux[ir].real() * chr->rho[is][ir] * sgn[ir*nspin+is];
+            }
+            vtxc += vtxc_vlapl;
+        }
     }
 
     //-------------------------------------------------
diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp
index b265af9bb15..ee66590b801 100644
--- a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp
+++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp
@@ -12,19 +12,17 @@
 //XC_POLARIZED, XC_UNPOLARIZED: internal flags used in LIBXC, denote the polarized(nspin=1) or unpolarized(nspin=2) calculations, definition can be found in xc.h from LIBXC
 void XC_Functional_Libxc::tau_xc(
             const std::vector<int> &func_id,
-    const double &rho, const double &grho, const double &atau, double &sxc,
-          double &v1xc, double &v2xc, double &v3xc)
+    const double &rho, const double &grho, const double &lapl_rho, const double &atau, double &sxc,
+          double &v1xc, double &v2xc, double &v3xc, double &vlaplc)
 {
-    double s, v1, v2, v3;
-	double lapl_rho, vlapl_rho;
-    lapl_rho = grho;
+    double s, v1, v2, v3, vlapl;
     std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED);
 
-    sxc = 0.0; v1xc = 0.0; v2xc = 0.0; v3xc = 0.0;
+    sxc = 0.0; v1xc = 0.0; v2xc = 0.0; v3xc = 0.0; vlaplc = 0.0;
 
     for(xc_func_type &func : funcs)
     {
-        xc_mgga_exc_vxc(&func,1,&rho,&grho,&lapl_rho,&atau,&s,&v1,&v2,&vlapl_rho,&v3);
+        xc_mgga_exc_vxc(&func,1,&rho,&grho,&lapl_rho,&atau,&s,&v1,&v2,&vlapl,&v3);
 #ifdef __EXX
         if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5)
         {
@@ -32,12 +30,14 @@ void XC_Functional_Libxc::tau_xc(
             v1 *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha);
             v2 *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha);
             v3 *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha);
+            vlapl *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha);
         }
 #endif
         sxc += s * rho;
         v2xc += v2 * 2.0;
         v1xc += v1;
         v3xc += v3;
+        vlaplc += vlapl;
     }
     XC_Functional_Libxc::finish_func(funcs);
 
@@ -49,16 +49,19 @@ void XC_Functional_Libxc::tau_xc_spin(
         const std::vector<int> &func_id,
         double rhoup, double rhodw,
         ModuleBase::Vector3<double> gdr1, ModuleBase::Vector3<double> gdr2,
+        double laplup, double lapldw,
         double tauup, double taudw,
         double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud,
-        double &v3xcup, double &v3xcdw)
+        double &v3xcup, double &v3xcdw, double &vlaplup, double &vlapldw)
 {
     sxc = v1xcup = v1xcdw = 0.0;
     v2xcup = v2xcdw = v2xcud = 0.0;
     v3xcup = v3xcdw = 0.0;
+    vlaplup = vlapldw = 0.0;
 
     const std::array<double,2> rho = {rhoup, rhodw};
     const std::array<double,3> grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()};
+    const std::array<double,2> lapl = {laplup, lapldw};
     const std::array<double,2> tau = {tauup, taudw};
 
 	std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED);
@@ -79,9 +82,9 @@ void XC_Functional_Libxc::tau_xc_spin(
             }
 
             double s = 0.0;
-            std::array<double,2> v1xc, v3xc, lapl, vlapl;
+            std::array<double,2> v1xc, v3xc, vlapl;
             std::array<double,3> v2xc;
-            // call Libxc function: xc_mgga_exc_vxc
+            // call Libxc function: xc_mgga_exc_vxc with real laplacian values
             xc_mgga_exc_vxc( &func, 1, rho.data(), grho.data(), lapl.data(), tau.data(), &s, v1xc.data(), v2xc.data(), vlapl.data(), v3xc.data());
 
             sxc += s * (rho[0] * sgn[0] + rho[1] * sgn[1]);
@@ -92,6 +95,8 @@ void XC_Functional_Libxc::tau_xc_spin(
             v2xcdw += 2.0 * v2xc[2] * sgn[1];
             v3xcup += v3xc[0] * sgn[0];
             v3xcdw += v3xc[1] * sgn[1];
+            vlaplup += vlapl[0] * sgn[0];
+            vlapldw += vlapl[1] * sgn[1];
         }
     }