diff --git a/examples/43-c6_coefficient.py b/examples/43-c6_coefficient.py
new file mode 100644
index 000000000..e2e05f427
--- /dev/null
+++ b/examples/43-c6_coefficient.py
@@ -0,0 +1,40 @@
+# Copyright 2021-2025 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+'''
+C6 coefficient
+'''
+
+import pyscf
+from gpu4pyscf import dft
+from gpu4pyscf.properties.c6 import calc_c6
+
+atom = '''
+O       0.0000000000    -0.0000000000     0.1174000000
+H      -0.7570000000    -0.0000000000    -0.4696000000
+H       0.7570000000     0.0000000000    -0.4696000000
+'''
+
+# Basis set must conatins diffuse functions and polarization functions
+bas='def2svpd'
+
+mol = pyscf.M(atom=atom, basis=bas)
+
+mf = dft.RKS(mol, xc='b3lyp').density_fit()
+e_gpu = mf.kernel() # -76.380311497689
+
+td_a = mf.TDDFT()
+td_b = mf.TDDFT()
+c6_val = calc_c6(td_a, td_b, n_grid=20) 
+print("Calculated C6 coefficient:", c6_val)
\ No newline at end of file
diff --git a/gpu4pyscf/properties/c6.py b/gpu4pyscf/properties/c6.py
new file mode 100644
index 000000000..baef5e2c2
--- /dev/null
+++ b/gpu4pyscf/properties/c6.py
@@ -0,0 +1,187 @@
+# Copyright 2021-2024 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import numpy as np
+import cupy as cp
+from pyscf import lib
+from gpu4pyscf.lib import logger
+from gpu4pyscf.lib.cupy_helper import contract
+from gpu4pyscf.tdscf import rhf as tdhf_gpu
+from gpu4pyscf.tdscf import rks as tdrks_gpu
+
+
+def _solve_full_spectrum(td):
+    log = logger.new_logger(td)
+    mf = td._scf
+    log.info('Constructing A and B matrices (GPU)...')
+    a_mat, b_mat = td.get_ab(mf)
+    
+    # Ensure matrices are real for RKS/RHF to avoid numerical noise creating complex components
+    a_mat = np.asarray(a_mat).real
+    b_mat = np.asarray(b_mat).real
+
+    nocc, nvir = a_mat.shape[:2]
+    nov = nocc * nvir
+
+    a_mat = a_mat.reshape(nov, nov)
+    b_mat = b_mat.reshape(nov, nov)
+    
+    is_tda = isinstance(td, (tdhf_gpu.TDA, tdrks_gpu.TDA))
+    
+    e_exc = None
+    xy_vectors = []
+
+    if is_tda:
+        log.info('Solving full TDA eigenvalue problem...')
+        w, x_mat = np.linalg.eigh(a_mat)
+        
+        for i in range(len(w)):
+            # TDA normalization
+            x_vec = x_mat[:, i].reshape(nocc, nvir)
+            x_vec *= np.sqrt(0.5)
+            xy_vectors.append((x_vec, 0))
+            
+        e_exc = w
+
+    else:
+        log.info('Solving full Casida eigenvalue problem...')
+        try:
+            amb = a_mat - b_mat
+            apb = a_mat + b_mat
+            
+            l_mat = np.linalg.cholesky(amb)
+            
+            # M = L.T * (A+B) * L
+            # M * v = w^2 * v
+            h_eff = np.dot(l_mat.T, np.dot(apb, l_mat))
+            w2, v = np.linalg.eigh(h_eff)
+            
+            mask = w2 > 1e-6
+            w_pos = np.sqrt(w2[mask])
+            v_pos = v[:, mask]
+            
+            # D = (L^T)^-1 * v
+            # Z = (1/w) * L * v 
+            d_vecs = np.linalg.solve(l_mat.T, v_pos)
+            z_vecs = np.dot(l_mat, v_pos) / w_pos[None, :]
+            
+            x_all = 0.5 * (z_vecs + d_vecs)
+            y_all = 0.5 * (z_vecs - d_vecs)
+            
+            for i in range(len(w_pos)):
+                x = x_all[:, i]
+                y = y_all[:, i]
+                
+                norm_x = np.linalg.norm(x)
+                norm_y = np.linalg.norm(y)
+                norm_diff = norm_x**2 - norm_y**2
+                
+                if abs(norm_diff) < 1e-9:
+                    scale = 1.0
+                else:
+                    scale = np.sqrt(0.5 / abs(norm_diff))
+                
+                x_vec = (x * scale).reshape(nocc, nvir)
+                y_vec = (y * scale).reshape(nocc, nvir)
+                
+                xy_vectors.append((x_vec, y_vec))
+            
+            e_exc = w_pos
+
+        except np.linalg.LinAlgError:
+            log.warn('Ground state unstable (A-B not positive definite). Fallback to non-symmetric diagonalization.')
+            
+            h_mat = np.empty((2 * nov, 2 * nov), dtype=a_mat.dtype)
+            h_mat[:nov, :nov] = a_mat
+            h_mat[:nov, nov:] = b_mat
+            h_mat[nov:, :nov] = -b_mat.conj()
+            h_mat[nov:, nov:] = -a_mat.conj()
+            
+            w, v = np.linalg.eig(h_mat)
+            
+            sorted_indices = np.argsort(w.real)
+            w = w[sorted_indices]
+            v = v[:, sorted_indices]
+            
+            mask = w.real > 1e-3
+            w_pos = w[mask]
+            v_pos = v[:, mask]
+            
+            for i in range(len(w_pos)):
+                xy_vec_c = v_pos[:, i]
+                idx_max = np.argmax(np.abs(xy_vec_c))
+                phase = np.angle(xy_vec_c[idx_max])
+                xy_vec = (xy_vec_c * np.exp(-1j * phase)).real
+                
+                x = xy_vec[:nov]
+                y = xy_vec[nov:]
+                
+                # Normalize: X^2 - Y^2 = 0.5 (PySCF convention for RHF/RKS)
+                norm_x = np.linalg.norm(x)
+                norm_y = np.linalg.norm(y)
+                norm_diff = norm_x**2 - norm_y**2
+                
+                if abs(norm_diff) < 1e-9:
+                    scale = 1.0
+                else:
+                    scale = np.sqrt(0.5 / abs(norm_diff))
+                
+                x_vec = (x * scale).reshape(nocc, nvir)
+                y_vec = (y * scale).reshape(nocc, nvir)
+                
+                xy_vectors.append((x_vec, y_vec))
+            
+            e_exc = w_pos.real
+
+    td.e = e_exc
+    td.xy = xy_vectors
+    td.converged = [True] * len(e_exc)
+    
+    return td
+
+def calc_c6(td_a, td_b, n_grid=20):
+    log = logger.new_logger(td_a)
+    log.info('\n' + '*' * 40)
+    log.info('GPU4PySCF C6 Calculation (Full Spectrum)')
+    log.info('*' * 40)
+    
+    x, w_leg = np.polynomial.legendre.leggauss(n_grid)
+    w0 = 0.5  # TODO: hard coded
+    freqs_im = w0 * (1 + x) / (1 - x)
+    weights = w_leg * w0 * 2 / ((1 - x)**2)
+    
+    log.info('Solving for System A')
+    _solve_full_spectrum(td_a)
+    f_osc_a = td_a.oscillator_strength()
+    e_exc_a = td_a.e
+    
+    log.info('Solving for System B')
+    _solve_full_spectrum(td_b)
+    f_osc_b = td_b.oscillator_strength() # in length gauge
+    e_exc_b = td_b.e
+    
+    # alpha(iw) = sum_I f_I / (w_I^2 + w^2)
+    denom_a = e_exc_a[:, None]**2 + freqs_im[None, :]**2
+    alpha_a = np.sum(f_osc_a[:, None] / denom_a, axis=0)
+    
+    denom_b = e_exc_b[:, None]**2 + freqs_im[None, :]**2
+    alpha_b = np.sum(f_osc_b[:, None] / denom_b, axis=0)
+    
+    integrand = alpha_a * alpha_b
+    c6_val = (3.0 / np.pi) * np.sum(integrand * weights)
+    
+    log.info(f'Calculated C6 coefficient: {c6_val:.6f} a.u.')
+    log.info('*' * 40 + '\n')
+    
+    return float(c6_val.real)
\ No newline at end of file
diff --git a/gpu4pyscf/properties/tests/test_c6.py b/gpu4pyscf/properties/tests/test_c6.py
new file mode 100644
index 000000000..c741251fb
--- /dev/null
+++ b/gpu4pyscf/properties/tests/test_c6.py
@@ -0,0 +1,171 @@
+# Copyright 2021-2025 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import unittest
+import numpy as np
+import cupy as cp
+from pyscf import lib, gto
+from gpu4pyscf import dft, tdscf
+from gpu4pyscf.properties import c6
+
+
+def diagonalize_casida(a, b, nroots=4):
+    nocc, nvir = a.shape[:2]
+    nov = nocc * nvir
+    a = a.reshape(nov, nov)
+    b = b.reshape(nov, nov)
+    h = np.block([[a        , b       ],
+                  [-b.conj(),-a.conj()]])
+    e = np.linalg.eig(np.asarray(h))[0]
+    lowest_e = np.sort(e[e.real > 0].real)[:nroots]
+    lowest_e = lowest_e[lowest_e > 1e-3]
+    return lowest_e
+
+
+def diagonalize_tda(a, nroots=5):
+    nocc, nvir = a.shape[:2]
+    nov = nocc * nvir
+    a = a.reshape(nov, nov)
+    e = np.linalg.eig(np.asarray(a))[0]
+    lowest_e = np.sort(e[e.real > 0].real)[:nroots]
+    lowest_e = lowest_e[lowest_e > 1e-3]
+    return lowest_e
+
+
+class KnownValues(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        mol = gto.Mole()
+        mol.verbose = 0
+        mol.output = '/dev/null'
+        mol.atom = [
+            ['H' , (0. , 0. , .917)],
+            ['F' , (0. , 0. , 0.)], ]
+        mol.basis = '631g'
+        cls.mol = mol.build()
+
+        cls.mf = cls.mol.RHF().to_gpu().run()
+        
+        mf_lda = cls.mol.RKS().to_gpu()
+        mf_lda.xc = 'lda, vwn'
+        cls.mf_lda = mf_lda.run(conv_tol=1e-10)
+
+        mf_b3lyp = cls.mol.RKS().to_gpu()
+        mf_b3lyp.xc = 'b3lyp'
+        cls.mf_b3lyp = mf_b3lyp.run(conv_tol=1e-10)
+
+        mol_h2o = gto.Mole()
+        mol_h2o.verbose = 0
+        mol_h2o.output = '/dev/null'
+        mol_h2o.atom = """
+        O 0.0000 0.0000 0.1173
+        H 0.0000 0.7572 -0.4692
+        H 0.0000 -0.7572 -0.4692
+        """
+        mol_h2o.basis = '631g'
+        cls.mol_h2o = mol_h2o.build()
+        cls.mf_h2o = cls.mol_h2o.RKS().to_gpu().run(xc='b3lyp')
+
+    @classmethod
+    def tearDownClass(cls):
+        cls.mol.stdout.close()
+        cls.mol_h2o.stdout.close()
+
+    def test_full_spectrum_tda_lda(self):
+        """Test full spectrum TDA solver against explicit diagonalization."""
+        td_b = self.mf_lda.TDA().set(nstates=5)
+        e_benchmark, xy_benchmark = td_b.kernel()
+        f_oscillator_benchmark = td_b.oscillator_strength()
+
+        td = self.mf_lda.TDA()
+        c6._solve_full_spectrum(td)
+        a, b = td.get_ab()
+        ref_e = diagonalize_tda(a, nroots=5)
+        f_oscillator = td.oscillator_strength()
+        
+        self.assertAlmostEqual(abs(td.e[:5] - ref_e).max(), 0, 6)
+        self.assertAlmostEqual(abs(td.e[:5] - e_benchmark).max(), 0, 6)
+        self.assertAlmostEqual(abs(f_oscillator[:5] - f_oscillator_benchmark).max(), 0, 6)
+        
+        x_vec = td.xy[0][0]
+        norm = 2.0 * np.sum(x_vec**2)
+        self.assertAlmostEqual(norm, 1.0, 5)
+
+    def test_full_spectrum_tddft_lda(self):
+        """Test full spectrum TDDFT solver (Pure) against explicit diagonalization."""
+        td_b = self.mf_lda.TDDFT().set(nstates=5)
+        e_benchmark, xy_benchmark = td_b.kernel()
+        f_oscillator_benchmark = td_b.oscillator_strength()
+
+        td = self.mf_lda.TDDFT()
+        c6._solve_full_spectrum(td)
+        f_oscillator = td.oscillator_strength()
+        
+        a, b = td.get_ab()
+        ref_e = diagonalize_casida(a, b, nroots=5)
+        
+        self.assertAlmostEqual(abs(td.e[:5] - ref_e).max(), 0, 6)
+        self.assertAlmostEqual(abs(td.e[:5] - e_benchmark).max(), 0, 6)
+        self.assertAlmostEqual(abs(f_oscillator[:5] - f_oscillator_benchmark).max(), 0, 6)
+
+    def test_full_spectrum_tddft_b3lyp(self):
+        """Test full spectrum TDDFT solver (Hybrid) against explicit diagonalization."""
+        td_b = self.mf_b3lyp.TDDFT().set(nstates=5)
+        e_benchmark, xy_benchmark = td_b.kernel()
+        f_oscillator_benchmark = td_b.oscillator_strength()
+
+        td = self.mf_b3lyp.TDDFT()
+        c6._solve_full_spectrum(td)
+        f_oscillator = td.oscillator_strength()
+        
+        a, b = td.get_ab()
+        ref_e = diagonalize_casida(a, b, nroots=5)
+        
+        self.assertAlmostEqual(abs(td.e[:5] - ref_e).max(), 0, 6)
+        self.assertAlmostEqual(abs(td.e[:5] - e_benchmark).max(), 0, 6)
+        self.assertAlmostEqual(abs(f_oscillator[:5] - f_oscillator_benchmark).max(), 0, 6)
+        
+        x_vec, y_vec = td.xy[0]
+        norm = np.sum(x_vec**2) - np.sum(y_vec**2)
+        self.assertAlmostEqual(norm, 0.5, 5)
+
+    def test_calc_c6_sanity(self):
+        """Test C6 calculation returns positive values and runs without error."""
+        td = self.mf_lda.TDDFT()
+        val = c6.calc_c6(td, td, n_grid=10)
+        ref_e = np.array([ 0.35545732,  0.35545732,  0.54368678,  1.1144098 ,  1.1144098 ,  1.1598993 ,
+            1.35171048,  1.41194799,  1.43205768,  1.43205769,  1.52953895,  1.52953895,
+            1.57505632,  1.57505632,  1.60406147,  1.90474313,  1.9660721 ,  1.9660721 ,
+            2.03275446,  2.09882741,  2.22934612,  2.22934612,  2.3315765 ,  2.91523625,
+            24.10618119, 24.86256929, 25.2256636 , 25.2256636 , 25.32903837, 25.77488066])
+
+        self.assertIsInstance(val, float)
+        self.assertAlmostEqual(abs(td.e - ref_e).max(), 0, 6)
+        self.assertAlmostEqual(val, 1.959496362180395, 6)
+        self.assertTrue(val > 0.0)
+
+    def test_calc_c6_symmetry(self):
+        """Test symmetry: C6(A, B) should equal C6(B, A)."""
+        td_a = self.mf_lda.TDA()     
+        td_b = self.mf_h2o.TDDFT()   
+        
+        c6_ab = c6.calc_c6(td_a, td_b, n_grid=12)
+        c6_ba = c6.calc_c6(td_b, td_a, n_grid=12)
+        
+        self.assertAlmostEqual(abs(c6_ab - c6_ba), 0, 9)
+        self.assertAlmostEqual(abs(c6_ab - 4.83381668540887), 0, 9)
+
+if __name__ == "__main__":
+    print("Full Tests for C6 calculations")
+    unittest.main()
\ No newline at end of file