diff --git a/source/Makefile.Objects b/source/Makefile.Objects
index ced2b16b823..d30e5198d44 100644
--- a/source/Makefile.Objects
+++ b/source/Makefile.Objects
@@ -866,6 +866,9 @@ OBJS_LCAO=evolve_elec.o\
       center2orb_orb22.o\
 
 OBJS_MODULE_RI=conv_coulomb_pot_k.o\
+      invmat_ri.o\
+      lri_cv.o\
+      lri_cv_tools.o\
       exx_abfs_jle.o\
       module_restart/restart_exx_csr.o\
       exx_abfs_io.o\
diff --git a/source/source_io/module_parameter/read_inp_exx_dftu.cpp b/source/source_io/module_parameter/read_inp_exx_dftu.cpp
index 2d60bdc7aaa..457afdd8bcf 100644
--- a/source/source_io/module_parameter/read_inp_exx_dftu.cpp
+++ b/source/source_io/module_parameter/read_inp_exx_dftu.cpp
@@ -264,7 +264,7 @@ void ReadInput::item_exx()
         item.annotation = "threshold to inverse Vq in abfs for generating Cs";
         item.category = "Exact Exchange (LCAO)";
         item.type = "Real";
-        item.description = "By default, the Coulomb matrix inversion required for obtaining LRI coefficients is performed using LU decomposition. However, this approach may suffer from numerical instabilities when a large set of auxiliary basis functions (ABFs) is employed. When exx_cs_inv_thr > 0, the inversion is instead carried out via matrix diagonalization. Eigenvalues smaller than exx_cs_inv_thr are discarded to improve numerical stability. A relatively safe and commonly recommended value is 1e-5.";
+        item.description = "By default, the Coulomb matrix inversion required for obtaining LRI coefficients is performed using Cholesky factorization. However, this approach may suffer from numerical instabilities when a large set of auxiliary basis functions (ABFs) is employed. When exx_cs_inv_thr > 0, the inversion is instead carried out via matrix diagonalization. Eigenvalues at or below exx_cs_inv_thr times the largest eigenvalue are discarded to improve numerical stability. A relatively safe and commonly recommended value is 1e-5. Force derivatives assume that modes remain on the same side of the cutoff. Eigenvalue crossings can make the truncated model nonsmooth during relaxation or molecular dynamics.";
         item.default_value = "-1";
         item.unit = "";
         read_sync_double(input.exx_cs_inv_thr);
diff --git a/source/source_lcao/module_ri/CMakeLists.txt b/source/source_lcao/module_ri/CMakeLists.txt
index d9373973845..9f1c903615a 100644
--- a/source/source_lcao/module_ri/CMakeLists.txt
+++ b/source/source_lcao/module_ri/CMakeLists.txt
@@ -5,6 +5,8 @@ if(ENABLE_LIBRI)
   list(
     APPEND
     objects
+    invmat_ri.cpp
+    lri_cv_tools.cpp
     matrix_orbs11.cpp
     matrix_orbs21.cpp
     matrix_orbs22.cpp
@@ -15,6 +17,7 @@ if(ENABLE_LIBRI)
     list(
       APPEND
       objects
+      lri_cv.cpp
       conv_coulomb_pot_k.cpp
       exx_abfs_ctor_orbs.cpp
       exx_abfs_io.cpp
diff --git a/source/source_lcao/module_ri/ewald_vq.hpp b/source/source_lcao/module_ri/ewald_vq.hpp
index d368b5565c7..d2c80047c68 100644
--- a/source/source_lcao/module_ri/ewald_vq.hpp
+++ b/source/source_lcao/module_ri/ewald_vq.hpp
@@ -6,6 +6,7 @@
 #ifndef EWALD_VQ_HPP
 #define EWALD_VQ_HPP
 
+#include "lri_cv_tools.h"
 #include <RI/comm/mix/Communicate_Tensors_Map_Judge.h>
 #include <RI/distribute/Distribute_Equally.h>
 #include <RI/global/Global_Func-1.h>
diff --git a/source/source_lcao/module_ri/exx_lri.hpp b/source/source_lcao/module_ri/exx_lri.hpp
index 618c58cef6f..83e2684aea2 100644
--- a/source/source_lcao/module_ri/exx_lri.hpp
+++ b/source/source_lcao/module_ri/exx_lri.hpp
@@ -6,6 +6,7 @@
 #ifndef EXX_LRI_HPP
 #define EXX_LRI_HPP
 
+#include "lri_cv_tools.h"
 #include "exx_lri.h"
 #include "source_cell/unitcell.h"
 #include "ri_2d_comm.h"
diff --git a/source/source_lcao/module_ri/invmat_ri.hpp b/source/source_lcao/module_ri/invmat_ri.cpp
similarity index 75%
rename from source/source_lcao/module_ri/invmat_ri.hpp
rename to source/source_lcao/module_ri/invmat_ri.cpp
index ebbe0208b65..4b83f1b1522 100644
--- a/source/source_lcao/module_ri/invmat_ri.hpp
+++ b/source/source_lcao/module_ri/invmat_ri.cpp
@@ -3,10 +3,9 @@
 // DATE :   2022-08-17
 //=======================
 
-#ifndef INVERSE_MATRIX_HPP
-#define INVERSE_MATRIX_HPP
 
 #include "invmat_ri.h"
+#include "abfs_construct_pca.h"
 #include "source_base/module_external/lapack_connector.h"
 
 #include <algorithm>
@@ -103,7 +102,7 @@ struct InverseMatrixTraits<std::complex<float>>
 };
 
 template <typename T>
-inline void Inverse_Matrix<T>::using_syev(const double& threshold_condition_number)
+void Inverse_Matrix<T>::using_syev(const double& threshold_condition_number)
 {
     using traits = InverseMatrixTraits<T>;
     using val_t = typename traits::value_type;
@@ -143,6 +142,71 @@ inline void Inverse_Matrix<T>::using_syev(const double& threshold_condition_numb
     BlasConnector::gemm(traits::gemm_trans, 'N', eA.nc, eA.nc, ie, 1, eA.c, eA.nc, eA.c, eA.nc, 0, A.ptr(), A.shape[1]);
 }
 
+template <typename T>
+std::array<RI::Tensor<T>, 3> Inverse_Matrix<T>::cal_inverse_with_derivatives(
+    const std::array<RI::Tensor<T>, 3>& derivatives,
+    const double threshold_condition_number)
+{
+    using traits = InverseMatrixTraits<T>;
+    using real_type = typename traits::value_type;
+    const std::size_t n = A.shape[0];
+    std::vector<real_type> eigenvalues(n);
+    int info = 0;
+    traits::syev(A, eigenvalues.data(), info);
+    if (info != 0)
+    {
+        throw std::range_error("Spectral inverse eigensolve failed, info=" + std::to_string(info));
+    }
+    const real_type maximum = *std::max_element(eigenvalues.begin(), eigenvalues.end());
+    const double cutoff = maximum * threshold_condition_number;
+    std::vector<bool> retained(n);
+    std::vector<real_type> inverse_values(n, 0);
+    for (std::size_t i = 0; i < n; ++i)
+    {
+        retained[i] = eigenvalues[i] > cutoff;
+        if (retained[i])
+        {
+            inverse_values[i] = 1 / eigenvalues[i];
+        }
+    }
+    // LAPACK eigenvectors are stored as rows U^H in the row-major tensor.
+    const RI::Tensor<T> vectors = A.copy();
+    const RI::Tensor<T> adjoint = vectors.conjugate().transpose();
+    RI::Tensor<T> weighted = vectors.copy();
+    for (std::size_t i = 0; i < n; ++i)
+    {
+        for (std::size_t j = 0; j < n; ++j)
+        {
+            weighted(i, j) *= inverse_values[i];
+        }
+    }
+    A = adjoint * weighted;
+    std::array<RI::Tensor<T>, 3> result;
+    for (std::size_t axis = 0; axis < 3; ++axis)
+    {
+        RI::Tensor<T> response = vectors * derivatives[axis] * adjoint;
+        for (std::size_t i = 0; i < n; ++i)
+        {
+            for (std::size_t j = 0; j < n; ++j)
+            {
+                real_type divided_difference = 0;
+                if (retained[i] && retained[j])
+                {
+                    divided_difference = -inverse_values[i] * inverse_values[j];
+                }
+                else if (retained[i] != retained[j])
+                {
+                    divided_difference = (inverse_values[i] - inverse_values[j])
+                                         / (eigenvalues[i] - eigenvalues[j]);
+                }
+                response(i, j) *= divided_difference;
+            }
+        }
+        result[axis] = adjoint * response * vectors;
+    }
+    return result;
+}
+
 /*
 void Inverse_Matrix::using_syev( const double &threshold_condition_number )
 {
@@ -259,4 +323,6 @@ void Inverse_Matrix<Tdata>::copy_down_triangle()
             A(i0, i1) = A(i1, i0);
 }
 
-#endif
+// Supported real and complex double-precision RI backends.
+template class Inverse_Matrix<double>;
+template class Inverse_Matrix<std::complex<double>>;
diff --git a/source/source_lcao/module_ri/invmat_ri.h b/source/source_lcao/module_ri/invmat_ri.h
index a1cb75cf1a9..6af02c40b06 100644
--- a/source/source_lcao/module_ri/invmat_ri.h
+++ b/source/source_lcao/module_ri/invmat_ri.h
@@ -7,10 +7,10 @@
 #define ABACUS_SOURCE_LCAO_MODULE_RI_INVMAT_RI_H
 
 
-#include "abfs_construct_pca.h"
 
 #include <RI/global/Tensor.h>
 #include <vector>
+#include <array>
 
 template <typename Tdata>
 class Inverse_Matrix
@@ -23,6 +23,11 @@ class Inverse_Matrix
     };
     void cal_inverse(const Method& method, const double& threshold_condition_number = 0.);
 
+    /// Invert with a fixed-rank spectral cutoff and differentiate the retained subspace.
+    std::array<RI::Tensor<Tdata>, 3> cal_inverse_with_derivatives(
+        const std::array<RI::Tensor<Tdata>, 3>& derivatives,
+        const double threshold_condition_number);
+
     void input(const RI::Tensor<Tdata>& m);
     void input(const std::vector<std::vector<RI::Tensor<Tdata>>>& ms);
     RI::Tensor<Tdata> output() const;
@@ -36,6 +41,5 @@ class Inverse_Matrix
     RI::Tensor<Tdata> A;
 };
 
-#include "invmat_ri.hpp"
 
 #endif // ABACUS_SOURCE_LCAO_MODULE_RI_INVMAT_RI_H
diff --git a/source/source_lcao/module_ri/lri_cv.hpp b/source/source_lcao/module_ri/lri_cv.cpp
similarity index 90%
rename from source/source_lcao/module_ri/lri_cv.hpp
rename to source/source_lcao/module_ri/lri_cv.cpp
index 49f197816c0..41d273d0b5b 100644
--- a/source/source_lcao/module_ri/lri_cv.hpp
+++ b/source/source_lcao/module_ri/lri_cv.cpp
@@ -3,8 +3,6 @@
 // DATE :   2022-08-17
 //=======================
 
-#ifndef LRI_CV_HPP
-#define LRI_CV_HPP
 
 #include "lri_cv.h"
 #include "lri_cv_tools.h"
@@ -406,11 +404,25 @@ LRI_CV<Tdata>::DPcal_C_dC(
                       DPcal_V(it1, it1, {0,0,0}, {{"writable_Vws",true}})}};
 
             std::vector<std::vector<RI::Tensor<Tdata>>> L;
+            std::array<LRI_CV_Tools::RIBlockMatrix<Tdata>, 3> dL;
             const double cs_inv_thr = this->p_info_ri != nullptr ? this->p_info_ri->Cs_inv_thr : 0.0;
-            if (cs_inv_thr > 0)
+            if (cs_inv_thr > 0 && !flag_finish_dC)
+            {
+                const auto derivative01 = DPcal_dV(it0, it1, R, flags);
+                const auto reverse_derivative = DPcal_dV(it1, it0, Rm, flags);
+                const auto derivative10 = LRI_CV_Tools::negative(reverse_derivative);
+                auto inverse = LRI_CV_Tools::cal_I_dI(V, derivative01, derivative10, cs_inv_thr);
+                L = std::move(inverse.first);
+                dL = std::move(inverse.second);
+            }
+            else if (cs_inv_thr > 0)
+            {
                 L = LRI_CV_Tools::cal_I(V, Inverse_Matrix<Tdata>::Method::syev, cs_inv_thr);
+            }
             else
+            {
                 L = LRI_CV_Tools::cal_I(V);
+            }
 
             const std::vector<RI::Tensor<Tdata>> C = LRI_CV_Tools::mul2(L,A);
             if(flags.at("writable_Cws"))
@@ -438,10 +450,26 @@ LRI_CV<Tdata>::DPcal_C_dC(
                                 this->index_abfs, this->index_lcaos, this->index_lcaos,
                                 Matrix_Orbs21::Matrix_Order::A1BA2))};
 
-                const std::array<RI::Tensor<Tdata>,3> dV_01 = DPcal_dV(it0, it1, R, flags);
-                const std::array<RI::Tensor<Tdata>,3> dV_10 = LRI_CV_Tools::negative(DPcal_dV(it1, it0, Rm, flags));
-
-                std::array<std::vector<RI::Tensor<Tdata>>,3>		// dC = L*(dA-dV*C)
+                std::array<std::vector<RI::Tensor<Tdata>>, 3> dC_tmp;
+                if (cs_inv_thr > 0)
+                {
+                    // C = F(V) A: include rotations between retained and discarded modes.
+                    for (std::size_t axis = 0; axis < 3; ++axis)
+                    {
+                        const std::vector<RI::Tensor<Tdata>> derivative_A = {dA[0][axis], dA[1][axis]};
+                        dC_tmp[axis] = LRI_CV_Tools::mul2(L, derivative_A);
+                        const auto inverse_response = LRI_CV_Tools::mul2(dL[axis], A);
+                        for (std::size_t block = 0; block < 2; ++block)
+                        {
+                            dC_tmp[axis][block] += inverse_response[block];
+                        }
+                    }
+                }
+                else
+                {
+                    const auto dV_01 = DPcal_dV(it0, it1, R, flags);
+                    const auto reverse_derivative = DPcal_dV(it1, it0, Rm, flags);
+                    const auto dV_10 = LRI_CV_Tools::negative(reverse_derivative);
                     dC_tmp = LRI_CV_Tools::mul2(
                             L,
                             LRI_CV_Tools::change_order( LRI_CV_Tools::minus(
@@ -449,6 +477,7 @@ LRI_CV<Tdata>::DPcal_C_dC(
                                 std::vector<std::array<RI::Tensor<Tdata>,3>>{
                                     LRI_CV_Tools::mul1(dV_01, C[1]),
                                     LRI_CV_Tools::mul1(dV_10, C[0])})));
+                }
                 const std::vector<std::array<RI::Tensor<Tdata>,3>>
                     dC = LRI_CV_Tools::change_order(std::move(dC_tmp));
                 if(flags.at("writable_dCws"))
@@ -464,5 +493,6 @@ LRI_CV<Tdata>::DPcal_C_dC(
     } // end else (!(C_read && flag_finish_dC))
 }
 
-
-#endif
+// Supported real and complex double-precision RI backends.
+template class LRI_CV<double>;
+template class LRI_CV<std::complex<double>>;
diff --git a/source/source_lcao/module_ri/lri_cv.h b/source/source_lcao/module_ri/lri_cv.h
index 12359d3b5e3..1103f6fcdd4 100644
--- a/source/source_lcao/module_ri/lri_cv.h
+++ b/source/source_lcao/module_ri/lri_cv.h
@@ -45,13 +45,13 @@ class LRI_CV
         const double &kmesh_times,
         std::shared_ptr<ORB_gaunt_table> MGT,
         const bool& init_C);
-    inline std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>
+    std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>
     cal_Vs(
         const UnitCell &ucell,
         const std::vector<TA> &list_A0,
         const std::vector<TAC> &list_A1,
         const std::map<std::string,bool> &flags);						// "writable_Vws"
-    inline std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>
+    std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>
     cal_dVs(
         const UnitCell &ucell,
         const std::vector<TA> &list_A0,
@@ -108,16 +108,16 @@ class LRI_CV
         const T_func_cal_Rcut& func_cal_Rcut,
         const T_func_DPcal_data<Tresult>& func_DPcal_data);
 
-    inline double cal_V_Rcut(const int it0, const int it1);
-    inline double cal_C_Rcut(const int it0, const int it1);
+    double cal_V_Rcut(const int it0, const int it1);
+    double cal_C_Rcut(const int it0, const int it1);
 
-    inline RI::Tensor<Tdata>
+    RI::Tensor<Tdata>
     DPcal_V(
         const int it0,
         const int it1,
         const Abfs::Vector3_Order<double> &R,
         const std::map<std::string,bool> &flags);						// "writable_Vws"
-    inline std::array<RI::Tensor<Tdata>,3>
+    std::array<RI::Tensor<Tdata>,3>
     DPcal_dV(
         const int it0,
         const int it1,
@@ -141,6 +141,5 @@ class LRI_CV
         const Tfunc &func_cal_o11);
 };
 
-#include "lri_cv.hpp"
 
 #endif
diff --git a/source/source_lcao/module_ri/lri_cv_tools.cpp b/source/source_lcao/module_ri/lri_cv_tools.cpp
new file mode 100644
index 00000000000..aadd628635f
--- /dev/null
+++ b/source/source_lcao/module_ri/lri_cv_tools.cpp
@@ -0,0 +1,51 @@
+#include "lri_cv_tools.h"
+
+#include <complex>
+
+template <typename T>
+std::pair<LRI_CV_Tools::RIBlockMatrix<T>, std::array<LRI_CV_Tools::RIBlockMatrix<T>, 3>>
+LRI_CV_Tools::cal_I_dI(const RIBlockMatrix<T>& metric,
+                      const std::array<RI::Tensor<T>, 3>& derivative01,
+                      const std::array<RI::Tensor<T>, 3>& derivative10,
+                      const double threshold)
+{
+    const std::vector<std::size_t> sizes = {metric[0][0].shape[0], metric[1][1].shape[0]};
+    std::array<RI::Tensor<T>, 3> derivatives;
+    for (std::size_t axis = 0; axis < 3; ++axis)
+    {
+        const RIBlockMatrix<T> blocks = {
+            {RI::Tensor<T>({sizes[0], sizes[0]}), derivative01[axis]},
+            {derivative10[axis], RI::Tensor<T>({sizes[1], sizes[1]})}};
+        Inverse_Matrix<T> pack;
+        pack.input(blocks);
+        derivatives[axis] = pack.output();
+    }
+    Inverse_Matrix<T> inverse;
+    inverse.input(metric);
+    const auto responses = inverse.cal_inverse_with_derivatives(derivatives, threshold);
+    std::array<RIBlockMatrix<T>, 3> block_responses;
+    for (std::size_t axis = 0; axis < 3; ++axis)
+    {
+        Inverse_Matrix<T> pack;
+        pack.input(responses[axis]);
+        block_responses[axis] = pack.output(sizes, sizes);
+    }
+    auto inverse_blocks = inverse.output(sizes, sizes);
+    return std::make_pair(inverse_blocks, block_responses);
+}
+
+template std::pair<LRI_CV_Tools::RIBlockMatrix<double>,
+                   std::array<LRI_CV_Tools::RIBlockMatrix<double>, 3>>
+LRI_CV_Tools::cal_I_dI(
+    const RIBlockMatrix<double>& metric,
+    const std::array<RI::Tensor<double>, 3>& derivative01,
+    const std::array<RI::Tensor<double>, 3>& derivative10,
+    const double threshold);
+
+template std::pair<LRI_CV_Tools::RIBlockMatrix<std::complex<double>>,
+                   std::array<LRI_CV_Tools::RIBlockMatrix<std::complex<double>>, 3>>
+LRI_CV_Tools::cal_I_dI(
+    const RIBlockMatrix<std::complex<double>>& metric,
+    const std::array<RI::Tensor<std::complex<double>>, 3>& derivative01,
+    const std::array<RI::Tensor<std::complex<double>>, 3>& derivative10,
+    const double threshold);
diff --git a/source/source_lcao/module_ri/lri_cv_tools.h b/source/source_lcao/module_ri/lri_cv_tools.h
index d5b0cd4f68f..9c9d2a48968 100644
--- a/source/source_lcao/module_ri/lri_cv_tools.h
+++ b/source/source_lcao/module_ri/lri_cv_tools.h
@@ -7,6 +7,7 @@
 #define LRI_CV_TOOLS_H
 
 #include "invmat_ri.h"
+#include "source_cell/unitcell.h"
 #include "abfs_vector3_order.h"
 #include "source_lcao/module_ri/abfs.h"
 
@@ -16,6 +17,7 @@
 #include <cstddef>
 #include <map>
 #include <vector>
+#include <utility>
 
 namespace LRI_CV_Tools
 {
@@ -30,6 +32,16 @@ extern std::vector<std::vector<RI::Tensor<Tdata>>> cal_I(const std::vector<std::
                                                          = Inverse_Matrix<Tdata>::Method::potrf,
                                                          const double& threshold_condition_number = 0.);
 
+template <typename T>
+using RIBlockMatrix = std::vector<std::vector<RI::Tensor<T>>>;
+
+template <typename T>
+std::pair<RIBlockMatrix<T>, std::array<RIBlockMatrix<T>, 3>> cal_I_dI(
+    const RIBlockMatrix<T>& metric,
+    const std::array<RI::Tensor<T>, 3>& derivative01,
+    const std::array<RI::Tensor<T>, 3>& derivative10,
+    const double threshold);
+
 template <typename Tdata>
 inline RI::Tensor<Tdata> transform_Rm(const RI::Tensor<Tdata>& V);
 template <typename Tdata>
diff --git a/source/source_lcao/module_ri/unittests/CMakeLists.txt b/source/source_lcao/module_ri/unittests/CMakeLists.txt
index e70aa0dbaf9..887eaa63c1b 100644
--- a/source/source_lcao/module_ri/unittests/CMakeLists.txt
+++ b/source/source_lcao/module_ri/unittests/CMakeLists.txt
@@ -13,7 +13,12 @@ if(ENABLE_LCAO)
     SOURCES
     test_exx_workspace.cpp
     ../exx_workspace.cpp
+    ../lri_cv.cpp
     ../../../source_basis/module_ao/orb_atomic_lm.cpp)
+  # Link the compiled LRI lifecycle while keeping unrelated integral assembly out
+  # of this focused geometry-sharing test.
+  target_compile_options(MODULE_RI_exx_workspace PRIVATE -ffunction-sections -fdata-sections)
+  target_link_options(MODULE_RI_exx_workspace PRIVATE -Wl,--gc-sections)
 endif()
 addtest(
   TARGET
@@ -37,3 +42,13 @@ addtest(
 addtest(TARGET MODULE_RI_abfs_vector3_order LIBS parameter SOURCES
         test_abfs_vector3_order.cpp)
 abacus_copy_test_files(support)
+
+addtest(
+  TARGET MODULE_RI_invmat_ri
+  LIBS parameter base device ${math_libs}
+  SOURCES test_invmat_ri.cpp ../invmat_ri.cpp ../lri_cv_tools.cpp
+          ../abfs_construct_pca.cpp
+          ../../../source_basis/module_ao/orb_atomic_lm.cpp)
+# Keep only the tensor eigensolver from the PCA translation unit in this focused test.
+target_compile_options(MODULE_RI_invmat_ri PRIVATE -ffunction-sections -fdata-sections)
+target_link_options(MODULE_RI_invmat_ri PRIVATE -Wl,--gc-sections)
diff --git a/source/source_lcao/module_ri/unittests/test_invmat_ri.cpp b/source/source_lcao/module_ri/unittests/test_invmat_ri.cpp
new file mode 100644
index 00000000000..0c31cad316c
--- /dev/null
+++ b/source/source_lcao/module_ri/unittests/test_invmat_ri.cpp
@@ -0,0 +1,244 @@
+#include "../invmat_ri.h"
+#include "../lri_cv_tools.h"
+#include "source_base/parallel_global.h"
+#include <gtest/gtest.h>
+#include <complex>
+
+namespace
+{
+template<typename T>
+void check_derivative(const double cutoff, const T phase)
+{
+    const double angle = 0.31;
+    const double step = 1e-5;
+    auto metric = [phase](const double x) {
+        RI::Tensor<T> q({3, 3});
+        q(0, 0) = std::cos(x);
+        q(0, 2) = -std::sin(x) * phase;
+        q(1, 1) = 1;
+        q(2, 0) = std::sin(x) * RI::Global_Func::get_conj(phase);
+        q(2, 2) = std::cos(x);
+        RI::Tensor<T> diagonal({3, 3});
+        diagonal(0, 0) = 2;
+        diagonal(1, 1) = 2;
+        diagonal(2, 2) = 0.01;
+        return q * diagonal * q.conjugate().transpose();
+    };
+    RI::Tensor<T> generator({3, 3});
+    generator(0, 2) = -phase;
+    generator(2, 0) = RI::Global_Func::get_conj(phase);
+    const auto v = metric(angle);
+    const auto dv = generator * v - v * generator;
+    const std::array<RI::Tensor<T>, 3> directions = {dv, RI::Tensor<T>({3, 3}), RI::Tensor<T>({3, 3})};
+    Inverse_Matrix<T> inverse;
+    inverse.input(v);
+    const auto derivatives = inverse.cal_inverse_with_derivatives(directions, cutoff);
+    Inverse_Matrix<T> plus;
+    const auto displaced_plus = metric(angle + step);
+    plus.input(displaced_plus);
+    plus.cal_inverse(Inverse_Matrix<T>::Method::syev, cutoff);
+    Inverse_Matrix<T> minus;
+    const auto displaced_minus = metric(angle - step);
+    minus.input(displaced_minus);
+    minus.cal_inverse(Inverse_Matrix<T>::Method::syev, cutoff);
+    const auto fd = (plus.output() - minus.output()) * T(0.5 / step);
+    const auto response = inverse.output() * v;
+    for (std::size_t i = 0; i < 9; ++i)
+    {
+        EXPECT_NEAR(std::abs(derivatives[0].ptr()[i] - fd.ptr()[i]), 0, 2e-6);
+        EXPECT_NEAR(std::abs(derivatives[1].ptr()[i]), 0, 1e-12);
+        EXPECT_NEAR(std::abs(derivatives[2].ptr()[i]), 0, 1e-12);
+    }
+    const double expected_trace = cutoff > 0 ? 2 : 3;
+    const T trace = response(0, 0) + response(1, 1) + response(2, 2);
+    EXPECT_NEAR(std::abs(trace - T(expected_trace)), 0, 1e-10);
+}
+}
+
+TEST(SpectralInverse, RealRetainedDiscardedRotation) { check_derivative<double>(0.1, 1.0); }
+TEST(SpectralInverse, ComplexRetainedDiscardedRotation)
+{
+    const auto phase = std::polar(1.0, 0.7);
+    check_derivative<std::complex<double>>(0.1, phase);
+}
+TEST(SpectralInverse, RealFullInverse) { check_derivative<double>(0.0, 1.0); }
+TEST(SpectralInverse, ComplexFullInverse)
+{
+    const auto phase = std::polar(1.0, 0.7);
+    check_derivative<std::complex<double>>(0.0, phase);
+}
+
+namespace
+{
+template<typename T>
+LRI_CV_Tools::RIBlockMatrix<T> two_center_metric(const std::array<double, 3>& x, const T phase)
+{
+    LRI_CV_Tools::RIBlockMatrix<T> metric = {
+        {RI::Tensor<T>({2, 2}), RI::Tensor<T>({2, 1})},
+        {RI::Tensor<T>({1, 2}), RI::Tensor<T>({1, 1})}};
+    metric[0][0](0, 0) = 3.0;
+    metric[0][0](0, 1) = 0.2;
+    metric[0][0](1, 0) = 0.2;
+    metric[0][0](1, 1) = 1.8;
+    metric[1][1](0, 0) = 0.8;
+    metric[0][1](0, 0) = 0.35 + 0.11 * x[0] - 0.07 * x[1] + 0.04 * x[2];
+    metric[0][1](1, 0) = phase * (0.18 - 0.06 * x[0] + 0.05 * x[1] + 0.09 * x[2]);
+    metric[1][0] = metric[0][1].conjugate().transpose();
+    return metric;
+}
+
+template<typename T>
+std::vector<RI::Tensor<T>> three_center_derivative(const std::size_t axis, const T phase)
+{
+    const std::array<double, 3> slope = {0.1, -0.08, 0.05};
+    std::vector<RI::Tensor<T>> derivative = {RI::Tensor<T>({2, 1, 2}), RI::Tensor<T>({1, 1, 2})};
+    for (std::size_t block = 0; block < 2; ++block)
+    {
+        for (std::size_t i = 0; i < derivative[block].shape[0]; ++i)
+        {
+            const std::size_t row = block == 0 ? i : 2;
+            for (std::size_t j = 0; j < 2; ++j)
+            {
+                derivative[block](i, 0, j) = phase * (slope[axis] * (row + 1) * (j + 0.5));
+            }
+        }
+    }
+    return derivative;
+}
+
+template<typename T>
+std::vector<RI::Tensor<T>> three_center(const std::array<double, 3>& x, const T phase)
+{
+    std::vector<RI::Tensor<T>> a = {RI::Tensor<T>({2, 1, 2}), RI::Tensor<T>({1, 1, 2})};
+    for (std::size_t block = 0; block < 2; ++block)
+    {
+        for (std::size_t i = 0; i < a[block].get_shape_all(); ++i)
+        {
+            a[block].ptr()[i] = T(0.3 + 0.1 * i + 0.2 * block);
+        }
+    }
+    for (std::size_t axis = 0; axis < 3; ++axis)
+    {
+        const auto derivative = three_center_derivative(axis, phase);
+        for (std::size_t block = 0; block < 2; ++block)
+        {
+            a[block] += derivative[block] * T(x[axis]);
+        }
+    }
+    return a;
+}
+
+template<typename T>
+void check_coefficients(const double cutoff, const T phase)
+{
+    const double step = 1e-5;
+    const std::array<double, 3> x = {0.13, -0.11, 0.07};
+    const auto metric = two_center_metric(x, phase);
+    const auto a = three_center(x, phase);
+    const std::array<double, 3> first_slope = {0.11, -0.07, 0.04};
+    const std::array<double, 3> second_slope = {-0.06, 0.05, 0.09};
+    std::array<RI::Tensor<T>, 3> derivative01;
+    std::array<RI::Tensor<T>, 3> derivative10;
+    for (std::size_t axis = 0; axis < 3; ++axis)
+    {
+        derivative01[axis] = RI::Tensor<T>({2, 1});
+        derivative01[axis](0, 0) = first_slope[axis];
+        derivative01[axis](1, 0) = phase * second_slope[axis];
+        derivative10[axis] = derivative01[axis].conjugate().transpose();
+    }
+    const auto inverse = LRI_CV_Tools::cal_I_dI(metric, derivative01, derivative10, cutoff);
+    const auto reference_inverse = LRI_CV_Tools::cal_I(metric, Inverse_Matrix<T>::Method::syev, cutoff);
+    for (std::size_t row = 0; row < 2; ++row)
+    {
+        for (std::size_t col = 0; col < 2; ++col)
+        {
+            EXPECT_EQ(inverse.first[row][col].shape[0], row == 0 ? 2 : 1);
+            EXPECT_EQ(inverse.first[row][col].shape[1], col == 0 ? 2 : 1);
+            for (std::size_t i = 0; i < inverse.first[row][col].get_shape_all(); ++i)
+            {
+                EXPECT_NEAR(std::abs(inverse.first[row][col].ptr()[i]
+                                     - reference_inverse[row][col].ptr()[i]), 0, 1e-12);
+            }
+        }
+    }
+    for (std::size_t axis = 0; axis < 3; ++axis)
+    {
+        auto plus_x = x;
+        auto minus_x = x;
+        plus_x[axis] += step;
+        minus_x[axis] -= step;
+        const auto plus_metric = two_center_metric(plus_x, phase);
+        const auto minus_metric = two_center_metric(minus_x, phase);
+        const auto plus_a = three_center(plus_x, phase);
+        const auto minus_a = three_center(minus_x, phase);
+        const auto plus_inverse = LRI_CV_Tools::cal_I(plus_metric, Inverse_Matrix<T>::Method::syev, cutoff);
+        const auto minus_inverse = LRI_CV_Tools::cal_I(minus_metric, Inverse_Matrix<T>::Method::syev, cutoff);
+        const auto plus_c = LRI_CV_Tools::mul2(plus_inverse, plus_a);
+        const auto minus_c = LRI_CV_Tools::mul2(minus_inverse, minus_a);
+        const auto da = three_center_derivative(axis, phase);
+        const auto first_term = LRI_CV_Tools::mul2(inverse.first, da);
+        const auto second_term = LRI_CV_Tools::mul2(inverse.second[axis], a);
+        for (std::size_t block = 0; block < 2; ++block)
+        {
+            const auto analytic = first_term[block] + second_term[block];
+            const auto fd = (plus_c[block] - minus_c[block]) * T(0.5 / step);
+            for (std::size_t i = 0; i < analytic.get_shape_all(); ++i)
+            {
+                EXPECT_NEAR(std::abs(analytic.ptr()[i] - fd.ptr()[i]), 0, 2e-8);
+            }
+        }
+        for (std::size_t row = 0; row < 2; ++row)
+        {
+            for (std::size_t col = 0; col < 2; ++col)
+            {
+                const auto fd = (plus_inverse[row][col] - minus_inverse[row][col]) * T(0.5 / step);
+                for (std::size_t i = 0; i < fd.get_shape_all(); ++i)
+                {
+                    EXPECT_NEAR(std::abs(inverse.second[axis][row][col].ptr()[i] - fd.ptr()[i]), 0, 2e-8);
+                }
+            }
+        }
+        // Fixed diagonal blocks and changing off-diagonal blocks vary the eigenvalues.
+        const double plus_norm = std::norm(plus_metric[0][1](0, 0)) + std::norm(plus_metric[0][1](1, 0));
+        const double minus_norm = std::norm(minus_metric[0][1](0, 0)) + std::norm(minus_metric[0][1](1, 0));
+        EXPECT_GT(std::abs(plus_norm - minus_norm), 1e-9);
+    }
+}
+}
+
+TEST(SpectralInverse, RealUnequalBlocksWithMovingThreeCenter)
+{
+    check_coefficients<double>(0.3, 1.0);
+}
+TEST(SpectralInverse, ComplexUnequalBlocksWithMovingThreeCenter)
+{
+    const auto phase = std::polar(1.0, 0.7);
+    check_coefficients<std::complex<double>>(0.3, phase);
+}
+TEST(SpectralInverse, RealUnequalBlocksWithoutTruncation)
+{
+    check_coefficients<double>(0.0, 1.0);
+}
+TEST(SpectralInverse, ComplexUnequalBlocksWithoutTruncation)
+{
+    const auto phase = std::polar(1.0, 0.7);
+    check_coefficients<std::complex<double>>(0.0, phase);
+}
+
+int main(int argc, char** argv)
+{
+    int processes = 1;
+    int threads = 1;
+    int rank = 0;
+    Parallel_Global::read_pal_param(argc, argv, processes, threads, rank);
+    POOL_WORLD = MPI_COMM_NULL;
+    KP_WORLD = MPI_COMM_NULL;
+    INT_BGROUP = MPI_COMM_NULL;
+    BP_WORLD = MPI_COMM_NULL;
+    GRID_WORLD = MPI_COMM_NULL;
+    DIAG_WORLD = MPI_COMM_NULL;
+    ::testing::InitGoogleTest(&argc, argv);
+    const int result = RUN_ALL_TESTS();
+    Parallel_Global::finalize_mpi();
+    return result;
+}