diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml
index 5aab944..2eb6daa 100644
--- a/.github/workflows/test.yaml
+++ b/.github/workflows/test.yaml
@@ -39,7 +39,7 @@ jobs:
         cd build;
         # run test on GPU node
         #srun -p GPU --time 1:0:0 --exclusive -N 1 -n 1 test/cu_solver_test
-        ctest -j 1
+        ctest -j 1 || ctest --rerun-failed --output-on-failure
 
     - name: Run coverage
       run: |
diff --git a/src/cu_symmetry.cu b/src/cu_symmetry.cu
index 980ab26..71ba46d 100644
--- a/src/cu_symmetry.cu
+++ b/src/cu_symmetry.cu
@@ -39,7 +39,7 @@ namespace green::gpu {
 
   cu_symmetry::~cu_symmetry() { release(); }
 
-  void cu_symmetry::initialize(const cu_symmetry_data& data, int nao, int naux, int nts, int ns) {
+  void cu_symmetry::initialize(const cu_symmetry_data& data, int nao, int nso, int naux, int nts, int ns) {
     release();
 
     const int nk  = static_cast<int>(data.nk);
@@ -48,11 +48,10 @@ namespace green::gpu {
     const int inq = static_cast<int>(data.inq);
 
     nao_  = nao;
+    nso_  = nso;
     naux_ = naux;
     nts_  = nts;
     ns_   = ns;
-    batch_count_   = static_cast<size_t>(nts) * static_cast<size_t>(ns);
-    matrix_stride_ = static_cast<size_t>(nao) * static_cast<size_t>(nao);
 
     // Cache host-side data for accessor methods and TR-conjugation check
     k_full_to_reduced_h_ = data.k_full_to_reduced;
@@ -81,26 +80,35 @@ namespace green::gpu {
     upload_array(k_tr_conj_d_,         data.k_tr_conj,         "k_tr_conj_d_");
     upload_array(q_tr_conj_d_,         data.q_tr_conj,         "q_tr_conj_d_");
 
+    // Determine matrix dimension of the AO-space transform: nao for scalar,
+    // nso (= 2·nao) for X2C. The same GEMM-based U·G·U† + TR-conj pipeline
+    // (transform_k_ao_device) handles both, with σ_y spinor mixing baked into
+    // the nso×nso U matrices coming from the input file.
     if (!data.k_ao_transforms.empty()) {
-      // Compute transform matrix dimension: nao for scalar, nso (= 2*nao) for X2C.
       size_t dim_sq = data.k_ao_transforms.size() / data.nk;
       k_transform_dim_ = static_cast<int>(std::round(std::sqrt(static_cast<double>(dim_sq))));
-      // Upload to GPU only for the scalar (nao×nao) case; X2C (nso×nso) uses transform_k_ao_device_2c instead.
-      if (k_transform_dim_ == nao_) {
-        std::vector<std::complex<float>> k_ao_f(data.k_ao_transforms.size());
-        cast_copy_complex(k_ao_f.data(), data.k_ao_transforms.data(), data.k_ao_transforms.size());
-
-        if (cudaMalloc(&k_ao_transform_full_d_, data.k_ao_transforms.size() * sizeof(cuDoubleComplex)) != cudaSuccess)
-          throw std::runtime_error("Failed to allocate k_ao_transform_full_d_.");
-        if (cudaMalloc(&k_ao_transform_full_f_, data.k_ao_transforms.size() * sizeof(cuComplex)) != cudaSuccess)
-          throw std::runtime_error("Failed to allocate k_ao_transform_full_f_.");
-        if (cudaMemcpy(k_ao_transform_full_d_, data.k_ao_transforms.data(), data.k_ao_transforms.size() * sizeof(std::complex<double>), cudaMemcpyHostToDevice) != cudaSuccess)
-          throw std::runtime_error("Failed to copy k_ao_transform_full to device.");
-        if (cudaMemcpy(k_ao_transform_full_f_, k_ao_f.data(), k_ao_f.size() * sizeof(std::complex<float>), cudaMemcpyHostToDevice) != cudaSuccess)
-          throw std::runtime_error("Failed to copy float k_ao_transform_full to device.");
-      }
+
+      std::vector<std::complex<float>> k_ao_f(data.k_ao_transforms.size());
+      cast_copy_complex(k_ao_f.data(), data.k_ao_transforms.data(), data.k_ao_transforms.size());
+
+      if (cudaMalloc(&k_ao_transform_full_d_, data.k_ao_transforms.size() * sizeof(cuDoubleComplex)) != cudaSuccess)
+        throw std::runtime_error("Failed to allocate k_ao_transform_full_d_.");
+      if (cudaMalloc(&k_ao_transform_full_f_, data.k_ao_transforms.size() * sizeof(cuComplex)) != cudaSuccess)
+        throw std::runtime_error("Failed to allocate k_ao_transform_full_f_.");
+      if (cudaMemcpy(k_ao_transform_full_d_, data.k_ao_transforms.data(), data.k_ao_transforms.size() * sizeof(std::complex<double>), cudaMemcpyHostToDevice) != cudaSuccess)
+        throw std::runtime_error("Failed to copy k_ao_transform_full to device.");
+      if (cudaMemcpy(k_ao_transform_full_f_, k_ao_f.data(), k_ao_f.size() * sizeof(std::complex<float>), cudaMemcpyHostToDevice) != cudaSuccess)
+        throw std::runtime_error("Failed to copy float k_ao_transform_full to device.");
+    } else {
+      k_transform_dim_ = nao;
     }
 
+    // Scratch sizing tracks the transform dimension: X2C input/output buffers
+    // are nso×nso per (k, t), so callers must pass ns=1 (X2C uses pseudo-ns=4
+    // for the 4 spin blocks but the matrix dim already absorbs spinor mixing).
+    batch_count_   = static_cast<size_t>(nts) * static_cast<size_t>(ns);
+    matrix_stride_ = static_cast<size_t>(k_transform_dim_) * static_cast<size_t>(k_transform_dim_);
+
     // q_p0_transforms stores U_q row-major (as-is). CUBLAS sees U_q^T (col-major).
     // Steps 2a/2c in compute_second_tau_contraction use OP_N/OP_C to recover U_q^T and U_q^*.
     if (!data.q_p0_transforms.empty()) {
@@ -190,12 +198,17 @@ namespace green::gpu {
     // Symmetry transform: G(k_full) = U * G(k_ibz) * U†  (non-TR)
     //                      G(k_full) = conj(U * G(k_ibz) * U†)  (TR)
     //
+    // Same pipeline handles scalar (dim = nao) and X2C (dim = nso = 2·nao): for X2C
+    // the σ_y spinor mixing is baked into the nso×nso U matrices read from the input
+    // file, so U·G·U† already produces correctly spin-flipped blocks at TR partners.
+    //
     // Both U and G are stored in ROW-MAJOR (ndarray / green-gpu MatrixXcd convention).
     // CUBLAS interprets row-major data as the transpose of the intended matrix.
     // To compute result_rm = U * G * U†, we need result_cm = (U*G*U†)^T = U^* * G^T * U^T.
     // With row-major U: OP_C → (U_rm^T)^H = U_rm^*, OP_N → U_rm^T.
     // With row-major G: OP_N → G_rm^T.
     // So: GEMM1(OP_C, OP_N) = U^* * G^T, GEMM2(OP_N, OP_N) = result * U^T.
+    const int dim = k_transform_dim_;
 
     if constexpr (std::is_same_v<cuda_complex_t, cuDoubleComplex>) {
       cuDoubleComplex* input_buf = input_scratch ? input_scratch : input_batch_z_d_;
@@ -209,15 +222,17 @@ namespace green::gpu {
       const cuDoubleComplex* U = k_ao_transform_full_d_ + k_full * matrix_stride_;
 
       // work = U^* * G^T  (col-major intermediate)
-      if (GEMM_STRIDED_BATCHED(handle, CUBLAS_OP_C, CUBLAS_OP_N, nao_, nao_, nao_, &one, U, nao_, 0, input_buf, nao_,
-                               matrix_stride_, &zero, work_buf, nao_, matrix_stride_, static_cast<int>(nts * ns)) != CUBLAS_STATUS_SUCCESS)
+      if (GEMM_STRIDED_BATCHED(handle, CUBLAS_OP_C, CUBLAS_OP_N, dim, dim, dim, &one, U, dim, 0, input_buf, dim,
+                               matrix_stride_, &zero, work_buf, dim, matrix_stride_, static_cast<int>(nts * ns)) != CUBLAS_STATUS_SUCCESS)
         throw std::runtime_error("Failed first batched GEMM in cu_symmetry::transform_k_ao_device_impl.");
       // out = work * U^T  (col-major result, row-major interpretation = U * G * U†)
-      if (GEMM_STRIDED_BATCHED(handle, CUBLAS_OP_N, CUBLAS_OP_N, nao_, nao_, nao_, &one, work_buf, nao_, matrix_stride_, U,
-                               nao_, 0, &zero, out_device, nao_, matrix_stride_, static_cast<int>(nts * ns)) != CUBLAS_STATUS_SUCCESS)
+      if (GEMM_STRIDED_BATCHED(handle, CUBLAS_OP_N, CUBLAS_OP_N, dim, dim, dim, &one, work_buf, dim, matrix_stride_, U,
+                               dim, 0, &zero, out_device, dim, matrix_stride_, static_cast<int>(nts * ns)) != CUBLAS_STATUS_SUCCESS)
         throw std::runtime_error("Failed second batched GEMM in cu_symmetry::transform_k_ao_device_impl.");
 
-      // TR: conjugate output to get conj(U * G_ibz * U†)
+      // TR: conjugate the full (U·G·U†) output to get conj(U * G_ibz * U†).
+      // For X2C this conjugates all four spin blocks together; the σ_y row/col
+      // mixing has already happened inside the GEMMs via the nso×nso U.
       if (k_tr_conj_h_.at(k_full) != 0) {
         const double minus_one = -1.0;
         if (RSCAL(handle, static_cast<int>(batch_elements), &minus_one, reinterpret_cast<double*>(out_device) + 1, 2) != CUBLAS_STATUS_SUCCESS)
@@ -235,15 +250,15 @@ namespace green::gpu {
       const cuComplex* U = k_ao_transform_full_f_ + k_full * matrix_stride_;
 
       // work = U^* * G^T  (col-major intermediate)
-      if (GEMM_STRIDED_BATCHED(handle, CUBLAS_OP_C, CUBLAS_OP_N, nao_, nao_, nao_, &one, U, nao_, 0, input_buf, nao_,
-                               matrix_stride_, &zero, work_buf, nao_, matrix_stride_, static_cast<int>(nts * ns)) != CUBLAS_STATUS_SUCCESS)
+      if (GEMM_STRIDED_BATCHED(handle, CUBLAS_OP_C, CUBLAS_OP_N, dim, dim, dim, &one, U, dim, 0, input_buf, dim,
+                               matrix_stride_, &zero, work_buf, dim, matrix_stride_, static_cast<int>(nts * ns)) != CUBLAS_STATUS_SUCCESS)
         throw std::runtime_error("Failed first batched GEMM in cu_symmetry::transform_k_ao_device_impl.");
       // out = work * U^T  (col-major result, row-major interpretation = U * G * U†)
-      if (GEMM_STRIDED_BATCHED(handle, CUBLAS_OP_N, CUBLAS_OP_N, nao_, nao_, nao_, &one, work_buf, nao_, matrix_stride_, U,
-                               nao_, 0, &zero, out_device, nao_, matrix_stride_, static_cast<int>(nts * ns)) != CUBLAS_STATUS_SUCCESS)
+      if (GEMM_STRIDED_BATCHED(handle, CUBLAS_OP_N, CUBLAS_OP_N, dim, dim, dim, &one, work_buf, dim, matrix_stride_, U,
+                               dim, 0, &zero, out_device, dim, matrix_stride_, static_cast<int>(nts * ns)) != CUBLAS_STATUS_SUCCESS)
         throw std::runtime_error("Failed second batched GEMM in cu_symmetry::transform_k_ao_device_impl.");
 
-      // TR: conjugate output to get conj(U * G_ibz * U†)
+      // TR: conjugate the full (U·G·U†) output to get conj(U * G_ibz * U†).
       if (k_tr_conj_h_.at(k_full) != 0) {
         const float minus_one = -1.0f;
         if (RSCAL(handle, static_cast<int>(batch_elements), &minus_one, reinterpret_cast<float*>(out_device) + 1, 2) != CUBLAS_STATUS_SUCCESS)
@@ -268,58 +283,4 @@ namespace green::gpu {
     transform_k_ao_device_impl(handle, stream, in_device, k_full, out_device, nts, ns, ibz_in_device, input_scratch, work_scratch);
   }
 
-  template <typename cuda_complex_t>
-  void cu_symmetry::transform_k_ao_device_2c_impl(cublasHandle_t handle, cudaStream_t stream,
-                                                  cuda_complex_t* ibz_in_device, size_t k_full,
-                                                  cuda_complex_t* out_device, int nts, int nao) {
-    using scalar_t = std::conditional_t<std::is_same_v<cuda_complex_t, cuDoubleComplex>, double, float>;
-    const size_t block_elems = static_cast<size_t>(nts) * nao * nao;
-    const size_t block_bytes = block_elems * sizeof(cuda_complex_t);
-
-    if (k_tr_conj_h_.at(k_full) == 0) {
-      // No TR: copy all 4 blocks unchanged.
-      cudaMemcpyAsync(out_device, ibz_in_device, 4 * block_bytes, cudaMemcpyDeviceToDevice, stream);
-    } else {
-      // TR needed (hardcoded minus_t=true block permutation):
-      //   ss=0 <- +conj(ibz ss=1)   [aa <- conj(bb)]
-      //   ss=1 <- +conj(ibz ss=0)   [bb <- conj(aa)]
-      //   ss=2 <- -conj(ibz ss=2)   [self]
-      //   ss=3 <- -conj(ibz ss=3)   [self]
-      // ibz_in_device and out_device are distinct, so ss=0/1 swap is safe.
-      cudaMemcpyAsync(out_device + 0 * block_elems, ibz_in_device + 1 * block_elems, block_bytes, cudaMemcpyDeviceToDevice, stream);
-      cudaMemcpyAsync(out_device + 1 * block_elems, ibz_in_device + 0 * block_elems, block_bytes, cudaMemcpyDeviceToDevice, stream);
-      cudaMemcpyAsync(out_device + 2 * block_elems, ibz_in_device + 2 * block_elems, block_bytes, cudaMemcpyDeviceToDevice, stream);
-      cudaMemcpyAsync(out_device + 3 * block_elems, ibz_in_device + 3 * block_elems, block_bytes, cudaMemcpyDeviceToDevice, stream);
-
-      cublasSetStream(handle, stream);
-      const scalar_t minus_one = static_cast<scalar_t>(-1.0);
-      // ss=0, ss=1: +conj -> negate imaginary parts only
-      if (RSCAL(handle, static_cast<int>(block_elems), &minus_one,
-                reinterpret_cast<scalar_t*>(out_device + 0 * block_elems) + 1, 2) != CUBLAS_STATUS_SUCCESS)
-        throw std::runtime_error("RSCAL conj ss=0 failed in transform_k_ao_device_2c.");
-      if (RSCAL(handle, static_cast<int>(block_elems), &minus_one,
-                reinterpret_cast<scalar_t*>(out_device + 1 * block_elems) + 1, 2) != CUBLAS_STATUS_SUCCESS)
-        throw std::runtime_error("RSCAL conj ss=1 failed in transform_k_ao_device_2c.");
-      // ss=2, ss=3: -conj -> negate real parts only
-      if (RSCAL(handle, static_cast<int>(block_elems), &minus_one,
-                reinterpret_cast<scalar_t*>(out_device + 2 * block_elems) + 0, 2) != CUBLAS_STATUS_SUCCESS)
-        throw std::runtime_error("RSCAL -conj ss=2 failed in transform_k_ao_device_2c.");
-      if (RSCAL(handle, static_cast<int>(block_elems), &minus_one,
-                reinterpret_cast<scalar_t*>(out_device + 3 * block_elems) + 0, 2) != CUBLAS_STATUS_SUCCESS)
-        throw std::runtime_error("RSCAL -conj ss=3 failed in transform_k_ao_device_2c.");
-    }
-  }
-
-  void cu_symmetry::transform_k_ao_device_2c(cublasHandle_t handle, cudaStream_t stream,
-                                              cuDoubleComplex* ibz_in_device, size_t k_full,
-                                              cuDoubleComplex* out_device, int nts, int nao) {
-    transform_k_ao_device_2c_impl(handle, stream, ibz_in_device, k_full, out_device, nts, nao);
-  }
-
-  void cu_symmetry::transform_k_ao_device_2c(cublasHandle_t handle, cudaStream_t stream,
-                                              cuComplex* ibz_in_device, size_t k_full,
-                                              cuComplex* out_device, int nts, int nao) {
-    transform_k_ao_device_2c_impl(handle, stream, ibz_in_device, k_full, out_device, nts, nao);
-  }
-
 }  // namespace green::gpu
\ No newline at end of file
diff --git a/src/cugw_qkpt.cu b/src/cugw_qkpt.cu
index e795c5c..df08d5b 100644
--- a/src/cugw_qkpt.cu
+++ b/src/cugw_qkpt.cu
@@ -325,10 +325,20 @@ namespace green::gpu {
           throw std::runtime_error("GEMM_STRIDED_BATCHED fails on gw_qkpt.compute_second_tau_contraction_2C().");
         }
         if (U_q != nullptr) {
+          // q-space symmetry transform: Y2 = U_q^left * P * U_q^right * Y1
           // U_q stored row-major as-is. CUBLAS sees U_q^T (col-major).
-          // q-space symmetry transform: U_q acts on auxiliary basis only, same as scalar case.
-          // 2a: effective op = OP_N(U_q^T) = U_q^T  →  T1 = U_q^T * Y1
-          if (GEMM_STRIDED_BATCHED(*handle_, CUBLAS_OP_N, CUBLAS_OP_T, naux_, nao2_, naux_, &one, U_q, naux_, 0,
+          //   OP_N(U_q^T) = U_q^T      OP_C(U_q^T) = U_q^*
+          // Non-TR: Y2 = U_q^T * P * U_q^* * Y1      (Left = OP_N, Right = OP_C)
+          // TR    : Y2 = U_q^* * P * U_q^T * Y1      (Left = OP_C, Right = OP_N)
+          // Folding the TR conjugation into the U_q OPs (mirroring the scalar
+          // path) is mathematically equivalent to applying conj(W) and avoids
+          // the post-step RSCAL on Y2, which would also conjugate Y1's
+          // contribution — Y1 already carries the correct TR convention from
+          // upstream (copy_Gk_2c on the CPU).
+          cublasOperation_t OP_Uq_Left  = q_conj_after_uq ? CUBLAS_OP_C : CUBLAS_OP_N;
+          cublasOperation_t OP_Uq_Right = q_conj_after_uq ? CUBLAS_OP_N : CUBLAS_OP_C;
+          // 2a: T1 = OP_Uq_Left(U_q^T) * Y1
+          if (GEMM_STRIDED_BATCHED(*handle_, OP_Uq_Left, CUBLAS_OP_T, naux_, nao2_, naux_, &one, U_q, naux_, 0,
                                    Y1t_Qin, nao2_, nauxnao2_, &zero, Y2t_inP, naux_, nauxnao2_,
                                    nt_mult) != CUBLAS_STATUS_SUCCESS) {
             throw std::runtime_error("GEMM_STRIDED_BATCHED fails on gw_qkpt.compute_second_tau_contraction_2C() [2a].");
@@ -339,21 +349,12 @@ namespace green::gpu {
                                    nt_mult) != CUBLAS_STATUS_SUCCESS) {
             throw std::runtime_error("GEMM_STRIDED_BATCHED fails on gw_qkpt.compute_second_tau_contraction_2C() [2b].");
           }
-          // 2c: effective op = OP_C(U_q^T) = U_q^*  →  Y2 = U_q^* * T2
-          if (GEMM_STRIDED_BATCHED(*handle_, CUBLAS_OP_C, CUBLAS_OP_N, naux_, nao2_, naux_, &one, U_q, naux_, 0,
+          // 2c: Y2 = OP_Uq_Right(U_q^T) * T2
+          if (GEMM_STRIDED_BATCHED(*handle_, OP_Uq_Right, CUBLAS_OP_N, naux_, nao2_, naux_, &one, U_q, naux_, 0,
                                    Y1t_Qin, naux_, nauxnao2_, &zero, Y2t_inP, naux_, nauxnao2_,
                                    nt_mult) != CUBLAS_STATUS_SUCCESS) {
             throw std::runtime_error("GEMM_STRIDED_BATCHED fails on gw_qkpt.compute_second_tau_contraction_2C() [2c].");
           }
-          // TR conjugation after U_q transform: conj(U * P * U†) = conjugate the result
-          if (q_conj_after_uq) {
-            scalar_t alpha = -1.0;
-            int      two   = 2;
-            if (RSCAL(*handle_, nauxnao2_ * nt_mult, &alpha, reinterpret_cast<scalar_t*>(Y2t_inP) + 1, two) !=
-                CUBLAS_STATUS_SUCCESS) {
-              throw std::runtime_error("RSCAL fails on gw_qkpt.compute_second_tau_contraction_2C() [q_conj_after_uq].");
-            }
-          }
         } else {
           // No q-space transform: Y2(Q,in) = P(Q,Q') * Y1^T(Q',in)
           if (GEMM_STRIDED_BATCHED(*handle_, CUBLAS_OP_N, CUBLAS_OP_T, naux_, nao2_, naux_, &one, Pqk_tQP + t * naux2_, naux_,
diff --git a/src/cugw_utils.cu b/src/cugw_utils.cu
index 0be1c36..b6a5fac 100644
--- a/src/cugw_utils.cu
+++ b/src/cugw_utils.cu
@@ -25,7 +25,7 @@ namespace green::gpu {
 
   template <typename prec>
   cugw_utils<prec>::cugw_utils(int _nts, int _nt_batch, int _nw_b, int _ns, int _nk, int _ink, int _nq, int _inq, int _nqkpt,
-                               int _NQ, int _nao, const cu_symmetry_data& sym_data, ztensor_view<5>& G_tskij_host,
+                               int _NQ, int _nao, int _nso, const cu_symmetry_data& sym_data, ztensor_view<5>& G_tskij_host,
                                bool low_device_memory, const MatrixXcd& Ttn_FB, const MatrixXcd& Tnt_BF,
                                LinearSolverType cuda_lin_solver, int _myid, int _intranode_rank, int _devCount_per_node) :
       _low_device_memory(low_device_memory), qkpts(_nqkpt), G_tskij_host_(G_tskij_host), V_Qpm(_NQ, _nao, _nao),
@@ -71,7 +71,7 @@ namespace green::gpu {
                                    g_kstij_device, g_ksmtij_device, sigma_kstij_device, sigma_k_locks);
     }
 
-    _cu_symmetry.initialize(sym_data, _nao, _NQ, _nts, _ns);
+    _cu_symmetry.initialize(sym_data, _nao, _nso, _NQ, _nts, _ns);
   }
 
   template <typename prec>
diff --git a/src/cuhf_utils.cu b/src/cuhf_utils.cu
index 113f39d..18d55a2 100644
--- a/src/cuhf_utils.cu
+++ b/src/cuhf_utils.cu
@@ -29,39 +29,51 @@ __global__ void initialize_array(cuDoubleComplex* array, cuDoubleComplex value,
 
 namespace green::gpu {
 
+  namespace {
+    // Common allocation shared by both constructors: V batch buffers, X/Y scratch,
+    // pinned host V buffer, and the per-k weights initialized to 1/nk.
+    void allocate_common_buffers(cuDoubleComplex** VkbatchQij, cuDoubleComplex** VkbatchaQj_conj,
+                                 cuDoubleComplex** X_kbatchQij, cuDoubleComplex** X_kbatchiaQ,
+                                 cuDoubleComplex** Y_kbatchij,  cuDoubleComplex** weights_fbz,
+                                 std::complex<double>** V_kQij_buffer,
+                                 cudaStream_t stream, size_t nk, size_t nkbatch, size_t NQnaosq, size_t naosq) {
+      using cuda_complex = cuDoubleComplex;
+      if (cudaMalloc(VkbatchQij, nkbatch * NQnaosq * sizeof(cuda_complex)) != cudaSuccess)
+        throw std::runtime_error("failure allocating Vkbatch");
+      if (cudaMalloc(VkbatchaQj_conj, nkbatch * NQnaosq * sizeof(cuda_complex)) != cudaSuccess)
+        throw std::runtime_error("failure allocating Vkbatch_conj");
+      if (cudaMallocHost(V_kQij_buffer, nkbatch * NQnaosq * sizeof(std::complex<double>)) != cudaSuccess)
+        throw std::runtime_error("failure allocating V on host");
+      if (cudaMalloc(X_kbatchQij, nkbatch * NQnaosq * sizeof(cuda_complex)) != cudaSuccess)
+        throw std::runtime_error("failure allocating XkbatchQij on device");
+      if (cudaMalloc(X_kbatchiaQ, nkbatch * NQnaosq * sizeof(cuda_complex)) != cudaSuccess)
+        throw std::runtime_error("failure allocating XkbatchiaQ on device");
+      if (cudaMalloc(Y_kbatchij, nkbatch * naosq * sizeof(cuda_complex)) != cudaSuccess)
+        throw std::runtime_error("failure allocating Y_kbatchij on device");
+      if (cudaMalloc(weights_fbz, nk * sizeof(cuda_complex)) != cudaSuccess)
+        throw std::runtime_error("failure allocating weights_fbz on device");
+      cuDoubleComplex nk_inv            = make_cuDoubleComplex(1. / nk, 0.);
+      int             threads_per_block = 512;
+      int             blocks_for_id     = nk / threads_per_block + 1;
+      initialize_array<<<blocks_for_id, threads_per_block, 0, stream>>>(*weights_fbz, nk_inv, nk);
+    }
+  }  // namespace
+
   cuhf_utils::cuhf_utils(size_t nk, size_t ink, size_t ns, size_t nao, size_t NQ, size_t nkbatch, ztensor<4> dm_fbz, int _myid,
                          int _intranode_rank, int _devCount_per_node) :
-      _nk(nk), _ink(ink), _ns(ns), _nao(nao), _NQ(NQ), _nkbatch(nkbatch), _naosq(nao * nao), _NQnaosq(NQ * nao * nao) {
+      _nao(nao), _nso(nao), _nsosq(0), _NQ(NQ), _naosq(nao * nao), _NQnaosq(NQ * nao * nao),
+      _ns(ns), _nk(nk), _ink(ink), _nkbatch(nkbatch) {
     if (cudaSetDevice(_intranode_rank % _devCount_per_node) != cudaSuccess) throw std::runtime_error("Error in cudaSetDevice1");
     if (cublasCreate(&_handle) != CUBLAS_STATUS_SUCCESS)
       throw std::runtime_error("Rank " + std::to_string(_myid) + ": error initializing cublas");
 
-    using cuda_complex = typename cu_type_map<std::complex<double>>::cuda_type;
     allocate_density_and_Fock(&_Dm_fbz_sk2ba, &_F_skij, dm_fbz.data(), _ink, _nk, _nao, _ns);
 
     if (cudaStreamCreate(&_stream) != cudaSuccess) throw std::runtime_error("main stream creation failed");
 
-    if (cudaMalloc(&_VkbatchQij, nkbatch * _NQnaosq * sizeof(cuda_complex)) != cudaSuccess)
-      throw std::runtime_error("failure allocating Vkbatch");
-    if (cudaMalloc(&_VkbatchaQj_conj, nkbatch * _NQnaosq * sizeof(cuda_complex)) != cudaSuccess)
-      throw std::runtime_error("failure allocating Vkbatch");
-
-    if (cudaMallocHost(&_V_kQij_buffer, nkbatch * _NQnaosq * sizeof(cxx_complex)) != cudaSuccess)
-      throw std::runtime_error("failure allocating V on host");
-
-    if (cudaMalloc(&_X_kbatchQij, _nkbatch * _NQnaosq * sizeof(cuda_complex)) != cudaSuccess)
-      throw std::runtime_error("failure allocating XkbatchQij on device");
-    if (cudaMalloc(&_X_kbatchiaQ, _nkbatch * _NQnaosq * sizeof(cuda_complex)) != cudaSuccess)
-      throw std::runtime_error("failure allocating XkbatchijQ on device");
-    if (cudaMalloc(&_Y_kbatchij, _nkbatch * _naosq * sizeof(cuda_complex)) != cudaSuccess)
-      throw std::runtime_error("failure allocating XkbatchijQ on device");
-
-    if (cudaMalloc(&_weights_fbz, _nk * sizeof(cuda_complex)) != cudaSuccess)
-      throw std::runtime_error("failure allocating weights_fbz on device");
-    cuDoubleComplex nk_inv            = make_cuDoubleComplex(1. / _nk, 0.);
-    int             threads_per_block = 512;
-    int             blocks_for_id     = _nk / threads_per_block + 1;
-    initialize_array<<<blocks_for_id, threads_per_block, 0, _stream>>>(_weights_fbz, nk_inv, _nk);
+    allocate_common_buffers(&_VkbatchQij, &_VkbatchaQj_conj, &_X_kbatchQij, &_X_kbatchiaQ, &_Y_kbatchij,
+                            &_weights_fbz, &_V_kQij_buffer, _stream, _nk, _nkbatch, _NQnaosq, _naosq);
+
     if (ns == 3) {
       _X2C = true;
     } else if (ns == 2 or ns == 1) {
@@ -71,6 +83,37 @@ namespace green::gpu {
     }
   }
 
+  cuhf_utils::cuhf_utils(size_t nk, size_t ink, size_t nao, size_t NQ, size_t nkbatch,
+                         const cuDoubleComplex* dm_fbz_nso_device,
+                         int _myid, int _intranode_rank, int _devCount_per_node) :
+      _X2C(true), _nao(nao), _nso(2 * nao), _nsosq((2 * nao) * (2 * nao)),
+      _NQ(NQ), _naosq(nao * nao), _NQnaosq(NQ * nao * nao),
+      _ns(3),  // pseudo-ns: aa, bb, ab spin blocks accumulated separately (ba derived later)
+      _nk(nk), _ink(ink), _nkbatch(nkbatch) {
+    if (cudaSetDevice(_intranode_rank % _devCount_per_node) != cudaSuccess) throw std::runtime_error("Error in cudaSetDevice1");
+    if (cublasCreate(&_handle) != CUBLAS_STATUS_SUCCESS)
+      throw std::runtime_error("Rank " + std::to_string(_myid) + ": error initializing cublas");
+
+    // X2C path keeps dm_fbz in device-resident (nk, nso, nso) layout. add_exchange_to_fock
+    // picks aa/bb/ab sub-views inside this buffer with lda=nso and per-ss offsets.
+    if (cudaMalloc(&_Dm_fbz_nso, _nk * _nsosq * sizeof(cuDoubleComplex)) != cudaSuccess)
+      throw std::runtime_error("failure allocating Dm_fbz_nso on device");
+    if (cudaMemcpy(_Dm_fbz_nso, dm_fbz_nso_device, _nk * _nsosq * sizeof(cuDoubleComplex),
+                   cudaMemcpyDeviceToDevice) != cudaSuccess)
+      throw std::runtime_error("failure copying Dm_fbz_nso input on device");
+
+    // F_skij still in 3-block (ss, ik, nao, nao) layout — downstream Fock assembly
+    // (copy_2c_Fock_from_device_to_host) expects this and derives ba = ab.adjoint().
+    if (cudaMalloc(&_F_skij, _ns * _ink * _naosq * sizeof(cuDoubleComplex)) != cudaSuccess)
+      throw std::runtime_error("failure allocating F_skij on device");
+    cudaMemset(_F_skij, 0, _ns * _ink * _naosq * sizeof(cuDoubleComplex));
+
+    if (cudaStreamCreate(&_stream) != cudaSuccess) throw std::runtime_error("main stream creation failed");
+
+    allocate_common_buffers(&_VkbatchQij, &_VkbatchaQj_conj, &_X_kbatchQij, &_X_kbatchiaQ, &_Y_kbatchij,
+                            &_weights_fbz, &_V_kQij_buffer, _stream, _nk, _nkbatch, _NQnaosq, _naosq);
+  }
+
   cuhf_utils::~cuhf_utils() {
     cudaStreamDestroy(_stream);
     cublasDestroy(_handle);
@@ -81,7 +124,8 @@ namespace green::gpu {
     cudaFree(_X_kbatchiaQ);
     cudaFree(_Y_kbatchij);
     cudaFree(_weights_fbz);
-    cudaFree(_Dm_fbz_sk2ba);
+    if (_Dm_fbz_sk2ba != nullptr) cudaFree(_Dm_fbz_sk2ba);
+    if (_Dm_fbz_nso   != nullptr) cudaFree(_Dm_fbz_nso);
     cudaFree(_F_skij);
 
     cudaFreeHost(_V_kQij_buffer);
@@ -113,10 +157,32 @@ namespace green::gpu {
     cublasSetStream(_handle, _stream);
 
     int nk_mult = std::min(_nkbatch, _nk - _k2);
+    const bool x2c_nso = (_Dm_fbz_nso != nullptr);
     for (size_t ss = 0; ss < _ns; ++ss) {
+      // Per-ss view into the dm_fbz buffer.
+      //   X2C nso layout: pick aa/bb/ab as nao×nao sub-views of the (k2, nso, nso)
+      //   block with lda=nso, stride=nsosq, plus a row/col offset selecting the spin
+      //   quadrant. ba (ss=3) is derived later as ab.adjoint().
+      //   Legacy 3-block (also non-X2C): contiguous (ss, k2, nao, nao), lda=nao,
+      //   stride=naosq.
+      cuda_complex* dm_ptr;
+      int           dm_lda;
+      long long     dm_stride;
+      if (x2c_nso) {
+        size_t row_off = (ss == 1) ? _nao : 0;       // ss=0 aa→0, ss=1 bb→nao, ss=2 ab→0
+        size_t col_off = (ss == 0) ? 0    : _nao;    // ss=0 aa→0, ss=1 bb→nao, ss=2 ab→nao
+        dm_ptr    = _Dm_fbz_nso + _k2 * _nsosq + row_off * _nso + col_off;
+        dm_lda    = static_cast<int>(_nso);
+        dm_stride = static_cast<long long>(_nsosq);
+      } else {
+        dm_ptr    = _Dm_fbz_sk2ba + ss * _nk * _naosq + _k2 * _naosq;
+        dm_lda    = static_cast<int>(_nao);
+        dm_stride = static_cast<long long>(_naosq);
+      }
+
       // X_skQia(k2, Qi, a) = VkbatchQij(k2, Qi, b) * Dm_fbz(s, k2, b, a)
       GEMM_STRIDED_BATCHED(_handle, CUBLAS_OP_N, CUBLAS_OP_N, _nao, _NQ * _nao, _nao, &one,
-                           _Dm_fbz_sk2ba + ss * _nk * _naosq + _k2 * _naosq, _nao, _naosq, _VkbatchQij, _nao, _NQnaosq, &zero,
+                           dm_ptr, dm_lda, dm_stride, _VkbatchQij, _nao, _NQnaosq, &zero,
                            _X_kbatchQij, _nao, _NQnaosq, nk_mult);
       // X_kbatchQij(k2, Q, ia) -> (k2, ia, Q)
       for (size_t kk2 = 0; kk2 < nk_mult; ++kk2) {
diff --git a/src/green/gpu/cu_symmetry.h b/src/green/gpu/cu_symmetry.h
index 7e04a60..ca89ca2 100644
--- a/src/green/gpu/cu_symmetry.h
+++ b/src/green/gpu/cu_symmetry.h
@@ -57,11 +57,11 @@ namespace green::gpu {
     std::vector<size_t> k1_from_k2q_map;   // k1 = k2 + q
     std::vector<size_t> k2_from_k1q_map;   // k2 = k1 - q
 
-    // Pre-built symmetry transforms (leave empty to skip GPU upload)
-    // For scalar: k_ao_transforms stores nao×nao matrices per k-point.
-    // For X2C: k_ao_transforms stores nso×nso (= 2nao × 2nao) matrices per k-point;
-    //          these are used for CPU-side transforms only (no GPU upload).
-    std::vector<std::complex<double>> k_ao_transforms;  // [nk][nao_or_nso][nao_or_nso]
+    // Pre-built symmetry transforms (leave empty to skip GPU upload).
+    // Scalar: nao×nao per k. X2C: nso×nso (= 2nao × 2nao) per k, with σ_y spinor
+    // mixing already encoded. Both are uploaded to the device and consumed by the
+    // same transform_k_ao_device GEMM pipeline.
+    std::vector<std::complex<double>> k_ao_transforms;  // [nk][dim][dim], dim ∈ {nao, nso}
     std::vector<std::complex<double>> q_p0_transforms;  // [nq][naux][naux]
   };
 
@@ -76,7 +76,7 @@ namespace green::gpu {
 
     ~cu_symmetry();
 
-    void initialize(const cu_symmetry_data& data, int nao, int naux, int nts, int ns);
+    void initialize(const cu_symmetry_data& data, int nao, int nso, int naux, int nts, int ns);
 
     bool initialized() const { return initialized_; }
 
@@ -98,14 +98,21 @@ namespace green::gpu {
       return q_p0_transform_full_f_ ? q_p0_transform_full_f_ + q_full * naux_ * naux_ : nullptr;
     }
 
-    // Transform device data already allocated (e.g., in qkpt buffers).
-    // Applies U_k * G * U_k^dagger on device; TR conjugation (conj(G)) is also applied on-device
-    // via RSCAL when the k-point is time-reversal related to its IBZ representative.
-    // Can optionally use a separate IBZ buffer as input (ibz_in_device) instead of in_device.
+    // Transform device data: applies U_k * G * U_k^dagger on contiguous matrices of size
+    // dim×dim, where dim is the AO-space dimension of the stored k_sym_transform_ao
+    // (nao for scalar runs, nso = 2·nao for X2C — the σ_y spinor mixing is baked into
+    // the nso×nso U). For TR-related k, conj(U * G * U†) is computed via a single
+    // RSCAL on the full output (no per-block manipulation).
     //
-    // input_scratch / work_scratch: optional per-caller scratch buffers, each nts*ns*nao*nao elements.
-    // When provided, these are used instead of the shared cu_symmetry scratch, enabling concurrent
-    // calls from different worker streams without data races.
+    // Layout expectations: in_device and out_device point to nts*ns matrices of size
+    // dim*dim each, contiguous and row-major (matching the green-gpu MatrixXcd
+    // convention). For X2C callers, ns should be 1 and the four spin blocks of G
+    // must be assembled into a single nso×nso block per (k, t) before the call.
+    //
+    // ibz_in_device: optional alternative source on device (used during chained ops).
+    // input_scratch / work_scratch: optional per-caller scratch buffers, each
+    // nts*ns*dim*dim elements. When provided, these are used instead of the shared
+    // cu_symmetry scratch — enabling concurrent calls from different worker streams.
     void transform_k_ao_device(cublasHandle_t handle, cudaStream_t stream, cuDoubleComplex* in_device, size_t k_full,
                                cuDoubleComplex* out_device, int nts, int ns,
                                cuDoubleComplex* ibz_in_device = nullptr,
@@ -117,24 +124,12 @@ namespace green::gpu {
                                cuComplex* input_scratch = nullptr,
                                cuComplex* work_scratch = nullptr);
 
-    // Device-side X2C TR spin-flip for G(k_ibz, -tau) → G(k_full, -tau).
-    // Input ibz_in_device holds the IBZ Green's function in 4-block layout [4, nts, nao, nao].
-    // If k_full is the IBZ representative (no TR needed), performs a direct device-to-device copy.
-    // If k_full is an anti-unitary image of its IBZ rep, applies the hardcoded minus_t=true
-    // block permutation + sign/conjugation: ss=0↔ss=1 with conj, ss=2,3 self with -conj.
-    // ibz_in_device and out_device must be distinct buffers.
-    void transform_k_ao_device_2c(cublasHandle_t handle, cudaStream_t stream,
-                                   cuDoubleComplex* ibz_in_device, size_t k_full,
-                                   cuDoubleComplex* out_device, int nts, int nao);
-    void transform_k_ao_device_2c(cublasHandle_t handle, cudaStream_t stream,
-                                   cuComplex* ibz_in_device, size_t k_full,
-                                   cuComplex* out_device, int nts, int nao);
-
   private:
     void release();
 
     bool initialized_ = false;
     int  nao_         = 0;
+    int  nso_         = 0;
     int  naux_        = 0;
     int  nts_         = 0;
     int  ns_          = 0;
@@ -146,11 +141,6 @@ namespace green::gpu {
                                     cuda_complex_t* out_device, int nts, int ns, cuda_complex_t* ibz_in_device,
                                     cuda_complex_t* input_scratch, cuda_complex_t* work_scratch);
 
-    template <typename cuda_complex_t>
-    void transform_k_ao_device_2c_impl(cublasHandle_t handle, cudaStream_t stream,
-                                       cuda_complex_t* ibz_in_device, size_t k_full,
-                                       cuda_complex_t* out_device, int nts, int nao);
-
     size_t*          k_full_to_reduced_d_   = nullptr;
     size_t*          k_reduced_to_full_d_   = nullptr;
     size_t*          q_full_to_reduced_d_   = nullptr;
diff --git a/src/green/gpu/cugw_utils.h b/src/green/gpu/cugw_utils.h
index 67bce30..16e8057 100644
--- a/src/green/gpu/cugw_utils.h
+++ b/src/green/gpu/cugw_utils.h
@@ -72,7 +72,7 @@ namespace green::gpu {
     using cuda_complex = typename cu_type_map<std::complex<prec>>::cuda_type;
 
   public:
-    cugw_utils(int _nts, int _nt_batch, int _nw_b, int _ns, int _nk, int _ink, int _nq, int _inq, int _nqkpt, int _NQ, int _nao,
+    cugw_utils(int _nts, int _nt_batch, int _nw_b, int _ns, int _nk, int _ink, int _nq, int _inq, int _nqkpt, int _NQ, int _nao, int _nso,
                const cu_symmetry_data& sym_data, ztensor_view<5>& G_tskij_host, bool _low_device_memory,
                const MatrixXcd& Ttn_FB, const MatrixXcd& Tnt_BF,
                LinearSolverType cuda_lin_solver, int _myid, int _intranode_rank, int _devCount_per_node);
diff --git a/src/green/gpu/cuhf_utils.h b/src/green/gpu/cuhf_utils.h
index e2456d4..77f8188 100644
--- a/src/green/gpu/cuhf_utils.h
+++ b/src/green/gpu/cuhf_utils.h
@@ -46,9 +46,19 @@ namespace green::gpu {
     using cuda_complex = typename cu_type_map<std::complex<double>>::cuda_type;
 
   public:
+    // Scalar (ns=1,2) and legacy X2C (ns=3, 3-block host dm_fbz layout) constructor.
     cuhf_utils(size_t nk, size_t ink, size_t ns, size_t nao, size_t NQ, size_t nkbatch, ztensor<4> dm_fbz, int _myid,
                int _intranode_rank, int _devCount_per_node);
 
+    // X2C constructor with device-resident dm_fbz in (nk, nso, nso) row-major layout.
+    // The per-spin-block GEMMs in add_exchange_to_fock pick aa/bb/ab sub-views with
+    // lda=nso and a per-ss (row, col) offset; ba is later derived as ab.adjoint() in
+    // copy_2c_Fock_from_device_to_host. Caller retains ownership of dm_fbz_nso_device;
+    // the data is copied into an internal device buffer at construction.
+    cuhf_utils(size_t nk, size_t ink, size_t nao, size_t NQ, size_t nkbatch,
+               const cuDoubleComplex* dm_fbz_nso_device,
+               int _myid, int _intranode_rank, int _devCount_per_node);
+
     ~cuhf_utils();
 
     static std::size_t size_divided_by_kbatch(size_t nao, size_t naux) {
@@ -78,6 +88,8 @@ namespace green::gpu {
 
     bool   _X2C;
     size_t _nao;
+    size_t _nso;     // 2*nao for X2C, == _nao for scalar (unused)
+    size_t _nsosq;   // _nso * _nso (X2C device-layout stride)
     size_t _NQ;
     size_t _naosq;
     size_t _NQnaosq;
@@ -88,7 +100,8 @@ namespace green::gpu {
     size_t _k2;
     size_t _k_pos;
     // Global objects
-    cuda_complex* _Dm_fbz_sk2ba;
+    cuda_complex* _Dm_fbz_sk2ba = nullptr;  // legacy 3-block (s, k2, b, a) layout
+    cuda_complex* _Dm_fbz_nso   = nullptr;  // X2C device-layout (k, nso, nso); non-null selects X2C path
     cuda_complex* _F_skij;
     cuda_complex* _weights_fbz;
     // Intermediate objects
diff --git a/src/green/gpu/gw_gpu_kernel.h b/src/green/gpu/gw_gpu_kernel.h
index 0edd8e8..f9fa4d8 100644
--- a/src/green/gpu/gw_gpu_kernel.h
+++ b/src/green/gpu/gw_gpu_kernel.h
@@ -32,10 +32,20 @@
 #include <array>
 
 #include "common_defs.h"
+#include "cu_symmetry.h"
 #include "df_integral_t.h"
 #include "gpu_kernel.h"
 
 namespace green::gpu {
+
+  // Build a cu_symmetry_data struct from a brillouin_zone_utils object.
+  // Called from g++-compiled code (HF or GW kernels); nvcc never sees HDF5 headers.
+  // build_k_ao: populate k_sym_transform_ao matrices (nao×nao for scalar, nso×nso for X2C).
+  // build_q_p0: populate q_sym_transform_p0 matrices (naux×naux). Set to false for HF.
+  cu_symmetry_data make_cu_symmetry_data(const symmetry::brillouin_zone_utils& bz,
+                                          int nao, int naux,
+                                          bool build_k_ao, bool build_q_p0);
+
   /**
    * @brief cuda GW Solver class performs self-energy calculation by means of GW approximation using density fitting
    */
diff --git a/src/gw_gpu_kernel.cpp b/src/gw_gpu_kernel.cpp
index 6f31bd1..5b83eb8 100644
--- a/src/gw_gpu_kernel.cpp
+++ b/src/gw_gpu_kernel.cpp
@@ -31,10 +31,11 @@
 namespace green::gpu {
   // Build a cu_symmetry_data struct from a brillouin_zone_utils object.
   // Called from g++-compiled code; never seen by nvcc.
-  static cu_symmetry_data make_cu_symmetry_data(const symmetry::brillouin_zone_utils& bz,
-                                                int nao, int naux,
-                                                bool build_k_ao, bool build_q_p0) {
+  cu_symmetry_data make_cu_symmetry_data(const symmetry::brillouin_zone_utils& bz,
+                                         int nao, int naux,
+                                         bool build_k_ao, bool build_q_p0) {
     cu_symmetry_data d;
+    const size_t nso  = bz.nso();
     const auto& ksym  = bz.k_symmetry();
     const auto& qsym  = bz.q_symmetry();
     const auto& kqmap = bz.k_q_map();
@@ -68,11 +69,11 @@ namespace green::gpu {
     }
 
     if (build_k_ao && nao > 0) {
-      d.k_ao_transforms.resize(d.nk * nao * nao);
-      MatrixXcd U_k(nao, nao);
+      d.k_ao_transforms.resize(d.nk * nso * nso);
+      MatrixXcd U_k(nso, nso);
       for (size_t k = 0; k < d.nk; ++k) {
         ksym.k_sym_transform_ao(U_k, k);
-        std::memcpy(d.k_ao_transforms.data() + k * nao * nao, U_k.data(), nao * nao * sizeof(std::complex<double>));
+        std::memcpy(d.k_ao_transforms.data() + k * nso * nso, U_k.data(), nso * nso * sizeof(std::complex<double>));
       }
     }
 
@@ -299,7 +300,7 @@ namespace green::gpu {
       // Build symmetry data on the CPU (g++ side) before passing to cugw_utils (nvcc side).
       // k-space AO transforms are only needed for scalar (non-relativistic) calculations.
       cu_symmetry_data sym_data = make_cu_symmetry_data(_bz_utils, _nao, _NQ, /*build_k_ao=*/true, /*build_q_p0=*/true);
-      cugw_utils<prec> cugw(_nts, _nt_batch, _nw_b, _ns, _nk, _ink, _nq, _inq, _nqkpt, _NQ, _nao, sym_data, g.object(),
+      cugw_utils<prec> cugw(_nts, _nt_batch, _nw_b, _ns, _nk, _ink, _nq, _inq, _nqkpt, _NQ, _nao, _nso, sym_data, g.object(),
                             _low_device_memory, _ft.Ttn_FB(), _ft.Tnt_BF(), _cuda_lin_solver, utils::context().global_rank,
                             utils::context().node_rank, _devCount_per_node);
       statistics.end();
@@ -468,9 +469,12 @@ namespace green::gpu {
       // Since the size of the Green's function and self-energy is 4 times largeer,
       // low_device_memory mode is always used.
       int pseudo_ns = 4;
-      // X2C: no k-space AO transforms needed; transform_k_ao_device_2c uses only TR flags.
-      cu_symmetry_data sym_data_x2c = make_cu_symmetry_data(_bz_utils, _nao, _NQ, /*build_k_ao=*/false, /*build_q_p0=*/true);
-      cugw_utils<prec> cugw(_nts, _nt_batch, _nw_b, pseudo_ns, _nk, _ink, _nq, _inq, _nqkpt, _NQ, _nao, sym_data_x2c,
+      // X2C: build nso×nso k-AO transforms with σ_y spinor mixing baked in.
+      // The same transform_k_ao_device pipeline that handles scalar U·G·U† + TR-conj
+      // also handles X2C — the only difference is the matrix dim (nso vs. nao),
+      // which cu_symmetry infers from the size of data.k_ao_transforms.
+      cu_symmetry_data sym_data_x2c = make_cu_symmetry_data(_bz_utils, _nao, _NQ, /*build_k_ao=*/true, /*build_q_p0=*/true);
+      cugw_utils<prec> cugw(_nts, _nt_batch, _nw_b, pseudo_ns, _nk, _ink, _nq, _inq, _nqkpt, _NQ, _nao, _nso, sym_data_x2c,
                             g.object(), true, _ft.Ttn_FB(), _ft.Tnt_BF(), _cuda_lin_solver,
                             utils::context().global_rank, utils::context().node_rank, _devCount_per_node);
       statistics.end();
@@ -526,54 +530,34 @@ namespace green::gpu {
       MPI_Win_unlock(0, sigma_tau.win());
     }
 
-    void x2c_gw_gpu_kernel::copy_Gk_2c(const ztensor<5> &G_tskij_host, tensor<std::complex<double>,4> &Gk_4tij, int k_full, bool minus_t) {
-      size_t k_ibz       = _bz_utils.k_symmetry().full_to_reduced()[k_full];
-      bool need_minus_k  = (_bz_utils.k_symmetry().tr_conj_list()[k_full] != 0);
-      for (size_t ss = 0; ss < 4; ++ss) {
-        size_t s1 = (ss % 2 == 0)? 0 : 1;
-        size_t s2 = ((ss+1) / 2 != 1)? 0 : 1;
-        size_t i_shift = (!minus_t)? s1*_nao : s2*_nao;
-        size_t j_shift = (!minus_t)? s2*_nao : s1*_nao;
-        for (size_t t = 0; t < _nts; ++t) {
-          size_t t_id = (!minus_t)? t : _nts-1-t;
-          if (!need_minus_k) {
-            matrix(Gk_4tij(ss, t)) = matrix(G_tskij_host(t_id, 0, k_ibz)).block(i_shift, j_shift, _nao, _nao);
-          } else {
-            size_t msi = (!minus_t)? (s1+1)%2 : (s2+1)%2;
-            size_t msj = (!minus_t)? (s2+1)%2 : (s1+1)%2;
-            if (msi == msj) {
-              matrix(Gk_4tij(ss, t)) = matrix(G_tskij_host(t_id, 0, k_ibz)).block(msi*_nao, msj*_nao, _nao, _nao).conjugate();
-            } else {
-              matrix(Gk_4tij(ss, t)) = (-1.0) * matrix(G_tskij_host(t_id, 0, k_ibz)).block(msi*_nao, msj*_nao, _nao, _nao).conjugate();
-            }
-          }
+    void x2c_gw_gpu_kernel::copy_Gk_2c(const ztensor<5>& G_tskij_host, tensor<std::complex<double>,4>& Gk_4tij, int k_full, bool minus_t) {
+      // Apply G(k) = U_k * G(ik) * U_k† via value_AO (handles both space group
+      // rotation and time-reversal conjugation via k_sym_transform_ao).
+      for (size_t t = 0; t < _nts; ++t) {
+        size_t    t_id = (!minus_t) ? t : _nts - 1 - t;
+        MatrixXcd G_k  = _bz_utils.k_symmetry().value_AO(G_tskij_host(t_id, 0), k_full);
+        for (size_t ss = 0; ss < 4; ++ss) {
+          size_t s1      = (ss % 2 == 0) ? 0 : 1;
+          size_t s2      = ((ss + 1) / 2 != 1) ? 0 : 1;
+          size_t i_shift = (!minus_t) ? s1 * _nao : s2 * _nao;
+          size_t j_shift = (!minus_t) ? s2 * _nao : s1 * _nao;
+          matrix(Gk_4tij(ss, t)) = G_k.block(i_shift, j_shift, _nao, _nao);
         }
       }
     }
 
-    void x2c_gw_gpu_kernel::copy_Gk_2c(const ztensor<5> &G_tskij_host, tensor<std::complex<float>,4> &Gk_4tij, int k_full, bool minus_t) {
-      size_t k_ibz       = _bz_utils.k_symmetry().full_to_reduced()[k_full];
-      bool need_minus_k  = (_bz_utils.k_symmetry().tr_conj_list()[k_full] != 0);
-      MatrixXcf G_ij(_nso, _nso);
-      for (size_t ss = 0; ss < 4; ++ss) {
-        size_t s1 = (ss % 2 == 0)? 0 : 1;
-        size_t s2 = ((ss+1) / 2 != 1)? 0 : 1;
-        size_t i_shift = (!minus_t)? s1*_nao : s2*_nao;
-        size_t j_shift = (!minus_t)? s2*_nao : s1*_nao;
-        for (size_t t = 0; t < _nts; ++t) {
-          size_t t_id = (!minus_t)? t : _nts-1-t;
-          G_ij = matrix(G_tskij_host(t_id, 0, k_ibz)).cast<std::complex<float> >();
-          if (!need_minus_k) {
-            matrix(Gk_4tij(ss, t)) = G_ij.block(i_shift, j_shift, _nao, _nao);
-          } else {
-            size_t msi = (!minus_t)? (s1+1)%2 : (s2+1)%2;
-            size_t msj = (!minus_t)? (s2+1)%2 : (s1+1)%2;
-            if (msi == msj) {
-              matrix(Gk_4tij(ss, t)) = G_ij.block(msi*_nao, msj*_nao, _nao, _nao).conjugate();
-            } else {
-              matrix(Gk_4tij(ss, t)) = (-1.0) * G_ij.block(msi*_nao, msj*_nao, _nao, _nao).conjugate();
-            }
-          }
+    void x2c_gw_gpu_kernel::copy_Gk_2c(const ztensor<5>& G_tskij_host, tensor<std::complex<float>,4>& Gk_4tij, int k_full, bool minus_t) {
+      // Apply G(k) = U_k * G(ik) * U_k† via value_AO (handles both space group
+      // rotation and time-reversal conjugation via k_sym_transform_ao).
+      for (size_t t = 0; t < _nts; ++t) {
+        size_t    t_id = (!minus_t) ? t : _nts - 1 - t;
+        MatrixXcf G_k  = _bz_utils.k_symmetry().value_AO(G_tskij_host(t_id, 0), k_full).cast<std::complex<float>>();
+        for (size_t ss = 0; ss < 4; ++ss) {
+          size_t s1      = (ss % 2 == 0) ? 0 : 1;
+          size_t s2      = ((ss + 1) / 2 != 1) ? 0 : 1;
+          size_t i_shift = (!minus_t) ? s1 * _nao : s2 * _nao;
+          size_t j_shift = (!minus_t) ? s2 * _nao : s1 * _nao;
+          matrix(Gk_4tij(ss, t)) = G_k.block(i_shift, j_shift, _nao, _nao);
         }
       }
     }
diff --git a/src/hf_gpu_kernel.cpp b/src/hf_gpu_kernel.cpp
index c9be09e..e196e0e 100644
--- a/src/hf_gpu_kernel.cpp
+++ b/src/hf_gpu_kernel.cpp
@@ -19,10 +19,18 @@
  * DEALINGS IN THE SOFTWARE.
  */
 
+#include <green/gpu/cu_symmetry.h>
 #include <green/gpu/cuhf_utils.h>
 #include <green/gpu/hf_gpu_kernel.h>
 
 namespace green::gpu {
+  // Defined in gw_gpu_kernel.cpp. Forward-declared here so we don't have to pull
+  // <green/gpu/gw_gpu_kernel.h> (with its grids/transformer dependencies) into the
+  // HF translation unit.
+  cu_symmetry_data make_cu_symmetry_data(const symmetry::brillouin_zone_utils& bz,
+                                         int nao, int naux,
+                                         bool build_k_ao, bool build_q_p0);
+
   void hf_gpu_kernel::HF_check_devices_free_space() {
     std::cout << std::setprecision(4) << std::boolalpha;
     // check devices' free space and determine nkbatch
@@ -227,15 +235,58 @@ namespace green::gpu {
 
   void x2c_hf_gpu_kernel::compute_exchange_selfenergy(ztensor<4> &new_Fock, const ztensor<4> &dm) {
     statistics.start("Initialization");
-    ztensor<4> dm_fbz_3kij(3, _nk, _nao, _nao);
-    get_dm_fbz(dm_fbz_3kij, dm);
-    // Also determines _nk_batch
+
+    // Build dm_fbz on device using the merged transform_k_ao_device pipeline:
+    //   dm_fbz(k_full) = U_k · dm_ibz(reduced_point(k_full)) · U_k†  (· conj if TR).
+    // The nso×nso U from the input file already encodes the σ_y spinor mixing for
+    // TR-related k, and the post-step RSCAL conjugates the full nso×nso result.
+    // Result lives on device in (nk, nso, nso) row-major layout; cuhf_utils picks
+    // per-spin-block (aa, bb, ab) sub-views with lda=nso at GEMM time.
+    cu_symmetry           sym;
+    cu_symmetry_data      sym_data = make_cu_symmetry_data(_bz_utils, _nao, /*naux=*/0,
+                                                            /*build_k_ao=*/true, /*build_q_p0=*/false);
+    sym.initialize(sym_data, _nao, _nso, /*naux=*/0, /*nts=*/1, /*ns=*/1);
+
+    const size_t     nsosq    = static_cast<size_t>(_nso) * _nso;
+    cuDoubleComplex* dm_ibz_d = nullptr;
+    cuDoubleComplex* dm_fbz_d = nullptr;
+    if (cudaMalloc(&dm_ibz_d, _ink * nsosq * sizeof(cuDoubleComplex)) != cudaSuccess)
+      throw std::runtime_error("Failed to allocate dm_ibz_d in x2c_hf_gpu_kernel::compute_exchange_selfenergy");
+    if (cudaMalloc(&dm_fbz_d, _nk * nsosq * sizeof(cuDoubleComplex)) != cudaSuccess)
+      throw std::runtime_error("Failed to allocate dm_fbz_d in x2c_hf_gpu_kernel::compute_exchange_selfenergy");
+
+    // dm has shape (ns=1, ink, nso, nso); the s=0 slice contiguously occupies the first ink*nsosq elements.
+    if (cudaMemcpy(dm_ibz_d, dm.data(), _ink * nsosq * sizeof(cuDoubleComplex), cudaMemcpyHostToDevice) != cudaSuccess)
+      throw std::runtime_error("Failed to upload dm_ibz to device");
+
+    cublasHandle_t local_handle;
+    if (cublasCreate(&local_handle) != CUBLAS_STATUS_SUCCESS)
+      throw std::runtime_error("Failed to create cublas handle for HF X2C transform");
+    cudaStream_t local_stream;
+    if (cudaStreamCreate(&local_stream) != cudaSuccess)
+      throw std::runtime_error("Failed to create cuda stream for HF X2C transform");
+
+    for (size_t k_full = 0; k_full < _nk; ++k_full) {
+      size_t k_ibz = sym.k_full_to_reduced(k_full);
+      sym.transform_k_ao_device(local_handle, local_stream,
+                                dm_ibz_d + k_ibz  * nsosq, k_full,
+                                dm_fbz_d + k_full * nsosq,
+                                /*nts=*/1, /*ns=*/1,
+                                /*ibz_in_device=*/nullptr,
+                                /*input_scratch=*/nullptr,
+                                /*work_scratch=*/nullptr);
+    }
+    cudaStreamSynchronize(local_stream);
+    cudaStreamDestroy(local_stream);
+    cublasDestroy(local_handle);
+    cudaFree(dm_ibz_d);
+
     HF_check_devices_free_space();
-    // Each process gets one cuda runner hf_utils
-    // Each NxN AO block of the 2-component exchange potential is evalulated individually
-    // using the non-relativistic functions with pseudo spin = 3 (i.e. aa, bb, ab blocks)
-    int pseudo_ns = 3;
-    cuhf_utils hf_utils(_nk, _ink, pseudo_ns, _nao, _NQ, _nk_batch, dm_fbz_3kij, utils::context().global_rank, utils::context().node_rank, _devCount_per_node);
+    // Each NxN AO block of the 2-component exchange potential is evaluated individually
+    // using pseudo-ns=3 (aa, bb, ab); ba is derived as ab.adjoint() in copy_2c_Fock_from_device_to_host.
+    cuhf_utils hf_utils(_nk, _ink, _nao, _NQ, _nk_batch, dm_fbz_d,
+                        utils::context().global_rank, utils::context().node_rank, _devCount_per_node);
+    cudaFree(dm_fbz_d);  // cuhf_utils owns its own copy now
     statistics.end();
     MPI_Barrier(_devices_comm);
 
@@ -302,61 +353,37 @@ namespace green::gpu {
 
   void x2c_hf_gpu_kernel::add_Ewald(ztensor<4>& new_Fock, const ztensor<4>& dm, const ztensor<4>& S, double madelung) {
     if (utils::context().global_rank < _ink * _ns) {
-      int direct_nprocs = (utils::context().global_size > _ink)? _ink : utils::context().global_size;
-      ztensor<3> dm_spblks[3] { {_ink, _nao, _nao}, {_ink, _nao, _nao}, {_ink, _nao, _nao} };
-      for (int ik = 0; ik < _ink; ++ik) {
-        CMMatrixXcd dmm(dm.data() + ik*_nso*_nso, _nso, _nso);
-        // alpha-alpha
-        matrix(dm_spblks[0](ik)) = dmm.block(0, 0, _nao, _nao);
-        // beta-beta
-        matrix(dm_spblks[1](ik)) = dmm.block(_nao, _nao, _nao, _nao);
-        // alpha-beta
-        matrix(dm_spblks[2](ik)) = dmm.block(0, _nao, _nao, _nao);
-      }
+      int direct_nprocs = (utils::context().global_size > _ink) ? _ink : utils::context().global_size;
+      // The Madelung uses the AO overlap S_aa (= S_bb) on both sides for all blocks —
+      // this is S_AO (spin-independent), not the spinor off-diagonal S_ab which is zero.
+      // s=0 (aa), s=1 (bb), s=2 (ab); ba is set as adjoint of the ab contribution.
       MatrixXcd buffer(_nao, _nao);
-      for (size_t iks = utils::context().global_rank; iks < 3*_ink; iks += direct_nprocs) {
-        size_t ik = iks / 3;
-        size_t is = iks % 3;
-        MMatrixXcd Fm_nso(new_Fock.data() + ik*_nso*_nso, _nso, _nso);
-        CMMatrixXcd Sm_nso(S.data() + ik*_nso*_nso, _nso, _nso);
-        MatrixXcd S_aa = Sm_nso.block(0, 0, _nao, _nao);
-        if (is == 0) {
-          // alpha-alpha
-          Fm_nso.block(0, 0, _nao, _nao) -= madelung * S_aa * matrix(dm_spblks[0](ik)) * S_aa;
-        } else if (is == 1) {
-          // beta-beta
-          Fm_nso.block(_nao, _nao, _nao, _nao) -= madelung * S_aa * matrix(dm_spblks[1](ik)) * S_aa;
-        } else if (is == 2) {
-          buffer = madelung * S_aa * matrix(dm_spblks[2](ik)) * S_aa;
-          // alpha-beta
-          Fm_nso.block(0, _nao, _nao, _nao) -= buffer;
-          // beta-alpha
+      for (size_t iks = utils::context().global_rank; iks < 3 * _ink; iks += direct_nprocs) {
+        size_t      ik  = iks / 3;
+        size_t      s   = iks % 3;
+        size_t      a   = s % 2;             // s=0→0 (alpha), s=1→1 (beta), s=2→0 (alpha)
+        size_t      b   = (s > 0) ? 1 : 0;  // s=0→0, s=1→1, s=2→1
+        MMatrixXcd  Fm_nso(new_Fock.data() + ik * _nso * _nso, _nso, _nso);
+        CMMatrixXcd Sm_nso(S.data() + ik * _nso * _nso, _nso, _nso);
+        CMMatrixXcd dm_nso(dm.data() + ik * _nso * _nso, _nso, _nso);
+        MatrixXcd   S_aa = Sm_nso.block(0, 0, _nao, _nao);  // AO overlap; same for all spin blocks
+        buffer = madelung * S_aa * dm_nso.block(a * _nao, b * _nao, _nao, _nao).eval() * S_aa;
+        Fm_nso.block(a * _nao, b * _nao, _nao, _nao) -= buffer;
+        if (a != b) {  // s=2 (ab): beta-alpha is adjoint of alpha-beta contribution
           Fm_nso.block(_nao, 0, _nao, _nao) -= buffer.conjugate().transpose();
         }
       }
     }
   }
 
-  void x2c_hf_gpu_kernel::get_dm_fbz(ztensor<4> &dm_fbz, const ztensor<4> &dm) {
-    size_t nsosq   = _nso*_nso;
+  void x2c_hf_gpu_kernel::get_dm_fbz(ztensor<4>& dm_fbz, const ztensor<4>& dm) {
+    // Apply G(k) = U_k * G(ik) * U_k† via value_AO (handles both space group
+    // rotation and time-reversal conjugation via k_sym_transform_ao).
     for (int k = 0; k < _nk; ++k) {
-      int         k_pos         = _bz_utils.k_symmetry().full_to_reduced()[k];
-      size_t shift_k = k_pos*nsosq;
-      CMMatrixXcd dmm(dm.data()+shift_k, _nso, _nso);
-
-      if (_bz_utils.k_symmetry().reduced_to_full()[k_pos] == k) {
-        // alpha-alpha
-        matrix(dm_fbz(0, k)) = dmm.block(0, 0, _nao, _nao);
-        // beta-beta
-        matrix(dm_fbz(1, k)) = dmm.block(_nao, _nao, _nao, _nao);
-        // alpha-beta
-        matrix(dm_fbz(2, k)) = dmm.block(0, _nao, _nao, _nao);
-      } else {
-        // -k counterpart (time-reversal Kramer pair)
-        matrix(dm_fbz(0, k)) = dmm.conjugate().block(_nao, _nao, _nao, _nao);
-        matrix(dm_fbz(1, k)) = dmm.conjugate().block(0, 0, _nao, _nao);
-        matrix(dm_fbz(2, k)) = (-1.0) * dmm.conjugate().block(_nao, 0, _nao, _nao);
-      }
+      MatrixXcd dm_k = _bz_utils.k_symmetry().value_AO(dm(0), k);
+      matrix(dm_fbz(0, k)) = dm_k.block(0,    0,    _nao, _nao);  // alpha-alpha
+      matrix(dm_fbz(1, k)) = dm_k.block(_nao, _nao, _nao, _nao);  // beta-beta
+      matrix(dm_fbz(2, k)) = dm_k.block(0,    _nao, _nao, _nao);  // alpha-beta
     }
   }
 }
diff --git a/test/cu_solver_test.cpp b/test/cu_solver_test.cpp
index f047f13..0e831e9 100644
--- a/test/cu_solver_test.cpp
+++ b/test/cu_solver_test.cpp
@@ -96,6 +96,163 @@ void solve_hf(const std::string& input, const std::string& int_hf, const std::st
 }
 
 
+// Check that one step of scf_type ("HF" or "GW") gives the same Sigma at every
+// IBZ k-point regardless of whether space-group, TR-only, or no symmetry is used.
+// Uses a cubic Ar (def2-svp, --x2c 2, 3x3x1 k-mesh) reference that exercises
+// both the orbital point-group representation (the basis includes p-shells)
+// and the SU(2) spinor part of the double-group transform under full X2C
+// spin-orbit coupling.  The 3x3x1 mesh gives non-trivial TR reduction
+// (no_symm=9 ink, trs_only=5 ink, full_symm=3 ink).  Tolerance accommodates
+// the integral-storage-floor disagreement between symmetry-reduced and
+// no-symmetry runs.
+void check_x2c_ar_symmetry(const std::string& scf_type, const std::string& lin, const std::string& mem) {
+  const std::string dir       = TEST_PATH + "/GW_X2C_Ar"s;
+  // Shared integral dir for both HF and GW (saves test data footprint).
+  const std::string df_path   = dir + "/df_hf_int"s;
+  const std::string grid_file = GRID_PATH + "/ir/1e4.h5"s;
+  constexpr double  tol = 1e-5;
+
+  size_t nts = 0;
+  {
+    green::h5pp::archive ar(grid_file);
+    ar["fermi/metadata/ncoeff"] >> nts;
+    nts += 2;
+  }
+
+  // Run one SCF step and return Sigma: [nts,ns,ink,nso,nso] for GW,
+  // [1,ns,ink,nso,nso] for HF (Sigma1 broadcast to uniform shape).
+  // All dimensions are read from input_file to avoid brittle hardcoding.
+  auto run = [&](const std::string& input_file, const std::string& data_file) {
+    size_t nao, nso, ns, nk, ink, NQ;
+    double madelung;
+    {
+      green::h5pp::archive ar(input_file);
+      ar["params/nao"]      >> nao;
+      ar["params/nso"]      >> nso;
+      ar["params/ns"]       >> ns;
+      ar["params/nk"]       >> nk;
+      ar["params/NQ"]       >> NQ;
+      ar["symmetry/k/ink"]  >> ink;
+      ar["HF/madelung"]     >> madelung;
+    }
+
+    auto        p    = green::params::params("DESCR");
+    std::string args = "test --restart 0 --itermax 1 --E_thr 1e-13 "
+                       "--mixing_type SIGMA_DAMPING --damping 0.7 "
+                       "--input_file=" + input_file + " --BETA 10 --grid_file=" + grid_file +
+                       " --cuda_low_gpu_memory " + mem + " --cuda_low_cpu_memory " + mem +
+                       (scf_type == "GW" ? " --dfintegral_file=" + df_path + " --cuda_linear_solver=" + lin
+                                         : " --dfintegral_hf_file=" + df_path);
+    green::grids::define_parameters(p);
+    green::gpu::custom_kernel_parameters(p);
+    green::symmetry::define_parameters(p);
+    p.define<std::string>("dfintegral_hf_file", "Path to HF integrals");
+    p.define<std::string>("dfintegral_file", "Path to integrals for correlated methods");
+    p.define<bool>("P_sp", "Compute polarization in single precision", false);
+    p.define<bool>("Sigma_sp", "Compute self-energy in single precision", false);
+    p.parse(args);
+    green::symmetry::brillouin_zone_utils bz(p);
+
+    green::gpu::ztensor<4> tmp(ns, nk, nso, nso), Sigma1(ns, ink, nso, nso), Sk(ns, ink, nso, nso);
+    {
+      green::h5pp::archive    ar(input_file);
+      green::gpu::dtensor<5>  S_k;
+      ar["HF/S-k"] >> S_k;
+      tmp << S_k.view<std::complex<double>>().reshape(ns, nk, nso, nso);
+    }
+    for (size_t is = 0; is < ns; ++is) Sk(is) << bz.full_to_ibz(tmp(is));
+
+    auto G_shared = green::utils::shared_object(green::gpu::ztensor<5>(nullptr, nts, ns, ink, nso, nso));
+    auto S_shared = green::utils::shared_object(green::gpu::ztensor<5>(nullptr, nts, ns, ink, nso, nso));
+    {
+      green::h5pp::archive ar(data_file, "r");
+      G_shared.fence();
+      if (!green::utils::context().node_rank) ar["G_tau"] >> G_shared.object();
+      G_shared.fence();
+    }
+
+    size_t                 nt_out = nts;
+    green::gpu::ztensor<5> result;
+    if (scf_type == "GW") {
+      green::grids::transformer_t ft(p);
+      auto [kernel, solver] = green::gpu::custom_gw_kernel(nso != nao, p, nao, nso, ns, NQ, ft, bz, Sk);
+      solver(G_shared, S_shared);
+      result.resize(nts, ns, ink, nso, nso);
+      S_shared.fence();
+      if (!green::utils::context().node_rank) result << S_shared.object();
+      S_shared.fence();
+    } else {
+      auto [kernel, solver] = green::gpu::custom_hf_kernel(nso != nao, p, nao, nso, ns, NQ, madelung, bz, Sk);
+      green::gpu::ztensor<4> dm(ns, ink, nso, nso);
+      dm << G_shared.object()(G_shared.object().shape()[0] - 1);
+      Sigma1 << solver(dm);
+      double prefactor = (ns == 2 or nao != nso) ? -1.0 : -2.0;
+      Sigma1 *= prefactor;
+      result.resize(1, ns, ink, nso, nso);
+      if (!green::utils::context().node_rank) result(0) << Sigma1;
+      nt_out = 1;
+    }
+    // Broadcast result to all ranks so assertions run on consistent data.
+    MPI_Bcast(result.data(), result.size(), MPI_CXX_DOUBLE_COMPLEX, 0, green::utils::context().global);
+    return std::make_pair(result, nt_out);
+  };
+
+  auto [Sigma_nosymm, nts_out] = run(dir + "/input_no_symm.h5"s,   dir + "/data_no_symm.h5"s);
+  auto [Sigma_symm,   _1]      = run(dir + "/input_full_symm.h5"s, dir + "/data_full_symm.h5"s);
+  auto [Sigma_trs,    _2]      = run(dir + "/input_trs_only.h5"s,  dir + "/data_trs_only.h5"s);
+
+  std::vector<long> ibz2bz_symm, ibz2bz_trs;
+  {
+    green::h5pp::archive ar(dir + "/input_full_symm.h5"s, "r");
+    ar["symmetry/k/ibz2bz"] >> ibz2bz_symm;
+  }
+  {
+    green::h5pp::archive ar(dir + "/input_trs_only.h5"s, "r");
+    ar["symmetry/k/ibz2bz"] >> ibz2bz_trs;
+  }
+
+  auto check = [&](const green::gpu::ztensor<5>& Sigma_ref, const green::gpu::ztensor<5>& Sigma_sym,
+                   const std::vector<long>& ibz2bz, const char* label) {
+    const size_t nso_ck = Sigma_sym.shape()[3];
+    const size_t nao_ck = nso_ck / 2;  // X2C: nso = 2 * nao
+    for (size_t i = 0; i < ibz2bz.size(); ++i) {
+      size_t k = ibz2bz[i];
+      for (size_t t = 0; t < nts_out; ++t) {
+        // Diagnostic: if any element of this (i,t) slice exceeds tol, dump
+        // per-spinor-block max diff before REQUIRE_THAT throws.
+        double overall_max = 0;
+        for (size_t r = 0; r < nso_ck; ++r)
+          for (size_t c = 0; c < nso_ck; ++c)
+            overall_max = std::max(overall_max,
+                std::abs(Sigma_sym(t, 0, i, r, c) - Sigma_ref(t, 0, k, r, c)));
+        if (overall_max > tol) {
+          auto block_max = [&](size_t r0, size_t c0, const char* name) {
+            double m = 0; size_t mr = 0, mc = 0;
+            for (size_t r = 0; r < nao_ck; ++r)
+              for (size_t c = 0; c < nao_ck; ++c) {
+                double d = std::abs(Sigma_sym(t, 0, i, r0 + r, c0 + c) -
+                                    Sigma_ref(t, 0, k, r0 + r, c0 + c));
+                if (d > m) { m = d; mr = r; mc = c; }
+              }
+            std::cout << "      " << name << " max=" << m
+                      << " at (" << mr << "," << mc << ")" << std::endl;
+          };
+          std::cout << "  [" << label << "] i=" << i << " k=" << k << " t=" << t
+                    << " overall max=" << overall_max << " (tol=" << tol << ")"
+                    << std::endl;
+          block_max(0,      0,      "aa");
+          block_max(nao_ck, nao_ck, "bb");
+          block_max(0,      nao_ck, "ab");
+          block_max(nao_ck, 0,      "ba");
+        }
+        REQUIRE_THAT(Sigma_sym(t, 0, i), IsCloseTo(Sigma_ref(t, 0, k), tol));
+      }
+    }
+  };
+  check(Sigma_nosymm, Sigma_symm, ibz2bz_symm, "full_symm vs no_symm");
+  check(Sigma_nosymm, Sigma_trs,  ibz2bz_trs,  "trs_only vs no_symm");
+}
+
 void solve_gw(const std::string& input, const std::string& int_f, const std::string& data, const std::string& lin, const std::string& mem, bool sp, const std::string& nt_batch) {
   auto        p           = green::params::params("DESCR");
   std::string input_file  = TEST_PATH + input;
@@ -210,6 +367,9 @@ TEST_CASE("GPU Solver") {
     solve_hf("/HF_X2C/input.h5", "/HF_X2C/df_hf_int", "/HF_X2C/data.h5", "true");
   }
 
+  SECTION("HF_X2C_Ar_Symmetry") { check_x2c_ar_symmetry("HF", "LU", "false"); }
+  SECTION("GW_X2C_Ar_Symmetry") { check_x2c_ar_symmetry("GW", "LU", "false"); }
+
   SECTION("Symmetry_Transform") {
     std::string input_file = TEST_PATH + "/GW/input.h5"s;
     auto p = green::params::params("DESCR");
diff --git a/test/cu_symmetry_test.cu b/test/cu_symmetry_test.cu
index 7e836b9..b8f2391 100644
--- a/test/cu_symmetry_test.cu
+++ b/test/cu_symmetry_test.cu
@@ -44,7 +44,7 @@ namespace green::gpu {
 
     // Initialize cu_symmetry on device
     cu_symmetry sym;
-    sym.initialize(sym_data, nao, /*naux=*/0, /*nts=*/1, /*ns=*/1);
+    sym.initialize(sym_data, nao, /*nso=*/nao, /*naux=*/0, /*nts=*/1, /*ns=*/1);
 
     // CUDA/CUBLAS setup
     cublasHandle_t handle;
diff --git a/test/data/GW_X2C_Ar/a.dat b/test/data/GW_X2C_Ar/a.dat
new file mode 100644
index 0000000..120a5a0
--- /dev/null
+++ b/test/data/GW_X2C_Ar/a.dat
@@ -0,0 +1,3 @@
+4.25    0.00    0.00
+0.00    4.25    0.00
+0.00    0.00    4.25
diff --git a/test/data/GW_X2C_Ar/atom.dat b/test/data/GW_X2C_Ar/atom.dat
new file mode 100644
index 0000000..1e2e4a0
--- /dev/null
+++ b/test/data/GW_X2C_Ar/atom.dat
@@ -0,0 +1 @@
+Ar  0.0 0.0 0.0
diff --git a/test/data/GW_X2C_Ar/data_full_symm.h5 b/test/data/GW_X2C_Ar/data_full_symm.h5
new file mode 100644
index 0000000..020640a
Binary files /dev/null and b/test/data/GW_X2C_Ar/data_full_symm.h5 differ
diff --git a/test/data/GW_X2C_Ar/data_no_symm.h5 b/test/data/GW_X2C_Ar/data_no_symm.h5
new file mode 100644
index 0000000..cc2bdbf
Binary files /dev/null and b/test/data/GW_X2C_Ar/data_no_symm.h5 differ
diff --git a/test/data/GW_X2C_Ar/data_trs_only.h5 b/test/data/GW_X2C_Ar/data_trs_only.h5
new file mode 100644
index 0000000..14817e9
Binary files /dev/null and b/test/data/GW_X2C_Ar/data_trs_only.h5 differ
diff --git a/test/data/GW_X2C_Ar/df_hf_int/VQ_0.h5 b/test/data/GW_X2C_Ar/df_hf_int/VQ_0.h5
new file mode 100644
index 0000000..700231c
Binary files /dev/null and b/test/data/GW_X2C_Ar/df_hf_int/VQ_0.h5 differ
diff --git a/test/data/GW_X2C_Ar/df_hf_int/meta.h5 b/test/data/GW_X2C_Ar/df_hf_int/meta.h5
new file mode 100644
index 0000000..89752d7
Binary files /dev/null and b/test/data/GW_X2C_Ar/df_hf_int/meta.h5 differ
diff --git a/test/data/GW_X2C_Ar/init.sh b/test/data/GW_X2C_Ar/init.sh
new file mode 100755
index 0000000..fe19990
--- /dev/null
+++ b/test/data/GW_X2C_Ar/init.sh
@@ -0,0 +1,26 @@
+#!/bin/bash
+
+
+export ScriptDir=$HOME/Documents/Codes/green-mbpt/python
+export basis="def2-svp"
+rm tmp*
+rm cderi*
+python $ScriptDir/init_data_df.py --x2c 2 --a a.dat --atom atom.dat \
+  --basis $basis --auxbasis def2-svp-jfit --nk 3 3 1 --xc PBE --keep_cderi true \
+  --space_symm true --tr_symm true --df_int 1 --output_path input_full_symm.h5 \
+  --use_j2c_eig_decomposition false
+
+rm tmp*
+rm cderi*
+python $ScriptDir/init_data_df.py --x2c 2 --a a.dat --atom atom.dat \
+  --basis $basis --auxbasis def2-svp-jfit --nk 3 3 1 --xc PBE --keep_cderi true \
+  --space_symm false --tr_symm true --df_int 0 --output_path input_trs_only.h5 \
+  --use_j2c_eig_decomposition false
+
+rm tmp*
+rm cderi*
+python $ScriptDir/init_data_df.py --x2c 2 --a a.dat --atom atom.dat \
+  --basis $basis --auxbasis def2-svp-jfit --nk 3 3 1 --xc PBE --keep_cderi true \
+  --space_symm false --tr_symm false --df_int 0 --output_path input_no_symm.h5 \
+  --use_j2c_eig_decomposition false
+
diff --git a/test/data/GW_X2C_Ar/input_full_symm.h5 b/test/data/GW_X2C_Ar/input_full_symm.h5
new file mode 100644
index 0000000..15e5309
Binary files /dev/null and b/test/data/GW_X2C_Ar/input_full_symm.h5 differ
diff --git a/test/data/GW_X2C_Ar/input_no_symm.h5 b/test/data/GW_X2C_Ar/input_no_symm.h5
new file mode 100644
index 0000000..06ab67c
Binary files /dev/null and b/test/data/GW_X2C_Ar/input_no_symm.h5 differ
diff --git a/test/data/GW_X2C_Ar/input_trs_only.h5 b/test/data/GW_X2C_Ar/input_trs_only.h5
new file mode 100644
index 0000000..c252a9b
Binary files /dev/null and b/test/data/GW_X2C_Ar/input_trs_only.h5 differ
diff --git a/test/data/GW_X2C_Ar/make_g_input.sh b/test/data/GW_X2C_Ar/make_g_input.sh
new file mode 100755
index 0000000..775f707
--- /dev/null
+++ b/test/data/GW_X2C_Ar/make_g_input.sh
@@ -0,0 +1,18 @@
+#!/bin/bash
+
+export GWBin=$HOME/Documents/Codes/green-mbpt/build/bin
+
+mpirun -n 3 $GWBin/mbpt.exe --BETA=10 --grid_file ir/1e4.h5 \
+  --input_file input_full_symm.h5 --results_file sim_full_symm.h5 \
+  --scf_type GW --mixing_weight 0.7 --itermax 1 --verbose 4 \
+  --dfintegral_file df_hf_int > out_gw_full_symm
+
+mpirun -n 3 $GWBin/mbpt.exe --BETA=10 --grid_file ir/1e4.h5 \
+  --input_file input_trs_only.h5 --results_file sim_trs_only.h5 \
+  --scf_type GW --mixing_weight 0.7 --itermax 1 --verbose 4 \
+  --dfintegral_file df_hf_int > out_gw_trs
+
+mpirun -n 3 $GWBin/mbpt.exe --BETA=10 --grid_file ir/1e4.h5 \
+  --input_file input_no_symm.h5 --results_file sim_no_symm.h5 \
+  --scf_type GW --mixing_weight 0.7 --itermax 1 --verbose 4 \
+  --dfintegral_file df_hf_int > out_gw_no_symm