From a384fb1d7be61a2b976e69ccd062d22cd8d787f2 Mon Sep 17 00:00:00 2001 From: "puzhichen.996" Date: Mon, 29 Dec 2025 16:40:33 +0800 Subject: [PATCH 1/7] add the c6 coefficients --- gpu4pyscf/properties/c6.py | 138 +++++++++++++++++++++ gpu4pyscf/properties/tests/test_c6.py | 165 ++++++++++++++++++++++++++ 2 files changed, 303 insertions(+) create mode 100644 gpu4pyscf/properties/c6.py create mode 100644 gpu4pyscf/properties/tests/test_c6.py diff --git a/gpu4pyscf/properties/c6.py b/gpu4pyscf/properties/c6.py new file mode 100644 index 000000000..d07d9be09 --- /dev/null +++ b/gpu4pyscf/properties/c6.py @@ -0,0 +1,138 @@ +# 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) + # a_mat = cp.asarray(a_mat) + # b_mat = cp.asarray(b_mat) + + 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) + + # w_cpu = cp.asnumpy(w) + # x_mat_cpu = cp.asnumpy(x_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...') + 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 + h_mat[nov:, nov:] = -a_mat + + 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 = v_pos[:, i] + 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(f'Solving for System A') + _solve_full_spectrum(td_a) + f_osc_a = td_a.oscillator_strength() + e_exc_a = td_a.e + + log.info(f'Solving for System B') + _solve_full_spectrum(td_b) + f_osc_b = td_b.oscillator_strength() + 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..1043103d6 --- /dev/null +++ b/gpu4pyscf/properties/tests/test_c6.py @@ -0,0 +1,165 @@ +# 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) + + self.assertIsInstance(val, float) + self.assertAlmostEqual(val, 1.9579916166911961, 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 From b9134562d9ab003c127fed803958e5c78d0d87ad Mon Sep 17 00:00:00 2001 From: "puzhichen.996" Date: Wed, 31 Dec 2025 09:52:45 +0800 Subject: [PATCH 2/7] add an example for c6 calculations. --- examples/43-c6_coefficient.py | 39 +++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 examples/43-c6_coefficient.py diff --git a/examples/43-c6_coefficient.py b/examples/43-c6_coefficient.py new file mode 100644 index 000000000..235b61892 --- /dev/null +++ b/examples/43-c6_coefficient.py @@ -0,0 +1,39 @@ +# 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 +''' + +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) # 45.027775533694815 \ No newline at end of file From e8346035a2c68471a0bad2ff4b01ad1244c4826e Mon Sep 17 00:00:00 2001 From: "puzhichen.996" Date: Wed, 31 Dec 2025 09:54:21 +0800 Subject: [PATCH 3/7] add some comments --- gpu4pyscf/properties/c6.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/gpu4pyscf/properties/c6.py b/gpu4pyscf/properties/c6.py index d07d9be09..096aef6b4 100644 --- a/gpu4pyscf/properties/c6.py +++ b/gpu4pyscf/properties/c6.py @@ -44,9 +44,6 @@ def _solve_full_spectrum(td): log.info('Solving full TDA eigenvalue problem...') w, x_mat = np.linalg.eigh(a_mat) - # w_cpu = cp.asnumpy(w) - # x_mat_cpu = cp.asnumpy(x_mat) - for i in range(len(w)): # TDA normalization x_vec = x_mat[:, i].reshape(nocc, nvir) @@ -119,7 +116,7 @@ def calc_c6(td_a, td_b, n_grid=20): log.info(f'Solving for System B') _solve_full_spectrum(td_b) - f_osc_b = td_b.oscillator_strength() + 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) From a700555a48eb52e80cd5bb786dd25ff8341d6787 Mon Sep 17 00:00:00 2001 From: "puzhichen.996" Date: Wed, 31 Dec 2025 11:20:20 +0800 Subject: [PATCH 4/7] fix a typo --- gpu4pyscf/properties/c6.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gpu4pyscf/properties/c6.py b/gpu4pyscf/properties/c6.py index 096aef6b4..8fb339e07 100644 --- a/gpu4pyscf/properties/c6.py +++ b/gpu4pyscf/properties/c6.py @@ -109,12 +109,12 @@ def calc_c6(td_a, td_b, n_grid=20): freqs_im = w0 * (1 + x) / (1 - x) weights = w_leg * w0 * 2 / ((1 - x)**2) - log.info(f'Solving for System A') + 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(f'Solving for System B') + 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 From 7d78fdf9fa9f43a97eae0e043bbc6bd4d01e07e9 Mon Sep 17 00:00:00 2001 From: "puzhichen.996" Date: Sun, 4 Jan 2026 09:18:39 +0800 Subject: [PATCH 5/7] add some debug codes --- gpu4pyscf/properties/tests/test_c6.py | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/gpu4pyscf/properties/tests/test_c6.py b/gpu4pyscf/properties/tests/test_c6.py index 1043103d6..59d3b5002 100644 --- a/gpu4pyscf/properties/tests/test_c6.py +++ b/gpu4pyscf/properties/tests/test_c6.py @@ -144,8 +144,28 @@ 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]) + ref_xy00 = np.array([[-1.30075599e-19+0.j, -2.46394827e-19+0.j, -3.57511290e-04+0.j, + 1.42099094e-04+0.j, -1.98997006e-18+0.j, -1.06832685e-18+0.j,], + [ 2.56736352e-17+0.j, -1.17104742e-16+0.j, -6.39696159e-03+0.j, + 2.54258390e-03+0.j, -1.13781911e-17+0.j, 1.70629105e-17+0.j,], + [ 9.53668504e-17+0.j, -1.18018049e-17+0.j, 8.06800482e-04+0.j, + -3.20676916e-04+0.j, -1.06770074e-16+0.j, -6.52977336e-17+0.j,], + [-1.42687004e-02+0.j, -2.33084878e-04+0.j, 1.95665208e-16+0.j, + 3.91504715e-17+0.j, 6.51825685e-05+0.j, -1.14309878e-04+0.j,], + [-7.06986911e-01+0.j, -1.15489115e-02+0.j, -1.07407259e-16+0.j, + -4.32672497e-17+0.j, 3.22967203e-03+0.j, -5.66383657e-03+0.j,],]) + + print("debug") + print(td.e) + print(td.xy[0][0]) self.assertIsInstance(val, float) + self.assertAlmostEqual(abs(td.e - ref_e).max(), 0, 6) + self.assertAlmostEqual(abs(td.xy[0][0] - ref_xy00).max(), 0, 6) self.assertAlmostEqual(val, 1.9579916166911961, 6) self.assertTrue(val > 0.0) From e9bc6056f1e7c65dd16d7ff7d112b5d10d187a42 Mon Sep 17 00:00:00 2001 From: "puzhichen.996" Date: Sun, 4 Jan 2026 14:35:16 +0800 Subject: [PATCH 6/7] fix a bug --- gpu4pyscf/properties/c6.py | 122 ++++++++++++++++++-------- gpu4pyscf/properties/tests/test_c6.py | 18 +--- 2 files changed, 89 insertions(+), 51 deletions(-) diff --git a/gpu4pyscf/properties/c6.py b/gpu4pyscf/properties/c6.py index 8fb339e07..baef5e2c2 100644 --- a/gpu4pyscf/properties/c6.py +++ b/gpu4pyscf/properties/c6.py @@ -26,8 +26,10 @@ def _solve_full_spectrum(td): mf = td._scf log.info('Constructing A and B matrices (GPU)...') a_mat, b_mat = td.get_ab(mf) - # a_mat = cp.asarray(a_mat) - # b_mat = cp.asarray(b_mat) + + # 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 @@ -54,43 +56,93 @@ def _solve_full_spectrum(td): else: log.info('Solving full Casida eigenvalue problem...') - 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 - h_mat[nov:, nov:] = -a_mat - - 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 = v_pos[:, i] - x = xy_vec[:nov] - y = xy_vec[nov:] + try: + amb = a_mat - b_mat + apb = a_mat + b_mat - # 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 + l_mat = np.linalg.cholesky(amb) - if abs(norm_diff) < 1e-9: - scale = 1.0 - else: - scale = np.sqrt(0.5 / abs(norm_diff)) + # 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) - x_vec = (x * scale).reshape(nocc, nvir) - y_vec = (y * scale).reshape(nocc, nvir) + mask = w2 > 1e-6 + w_pos = np.sqrt(w2[mask]) + v_pos = v[:, mask] - xy_vectors.append((x_vec, y_vec)) - - e_exc = w_pos.real + # 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 diff --git a/gpu4pyscf/properties/tests/test_c6.py b/gpu4pyscf/properties/tests/test_c6.py index 59d3b5002..c741251fb 100644 --- a/gpu4pyscf/properties/tests/test_c6.py +++ b/gpu4pyscf/properties/tests/test_c6.py @@ -149,24 +149,10 @@ def test_calc_c6_sanity(self): 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]) - ref_xy00 = np.array([[-1.30075599e-19+0.j, -2.46394827e-19+0.j, -3.57511290e-04+0.j, - 1.42099094e-04+0.j, -1.98997006e-18+0.j, -1.06832685e-18+0.j,], - [ 2.56736352e-17+0.j, -1.17104742e-16+0.j, -6.39696159e-03+0.j, - 2.54258390e-03+0.j, -1.13781911e-17+0.j, 1.70629105e-17+0.j,], - [ 9.53668504e-17+0.j, -1.18018049e-17+0.j, 8.06800482e-04+0.j, - -3.20676916e-04+0.j, -1.06770074e-16+0.j, -6.52977336e-17+0.j,], - [-1.42687004e-02+0.j, -2.33084878e-04+0.j, 1.95665208e-16+0.j, - 3.91504715e-17+0.j, 6.51825685e-05+0.j, -1.14309878e-04+0.j,], - [-7.06986911e-01+0.j, -1.15489115e-02+0.j, -1.07407259e-16+0.j, - -4.32672497e-17+0.j, 3.22967203e-03+0.j, -5.66383657e-03+0.j,],]) - - print("debug") - print(td.e) - print(td.xy[0][0]) + self.assertIsInstance(val, float) self.assertAlmostEqual(abs(td.e - ref_e).max(), 0, 6) - self.assertAlmostEqual(abs(td.xy[0][0] - ref_xy00).max(), 0, 6) - self.assertAlmostEqual(val, 1.9579916166911961, 6) + self.assertAlmostEqual(val, 1.959496362180395, 6) self.assertTrue(val > 0.0) def test_calc_c6_symmetry(self): From 1596684bf11f4865b310252e4d3ef1eb54bf4b47 Mon Sep 17 00:00:00 2001 From: "puzhichen.996" Date: Tue, 6 Jan 2026 12:39:07 +0800 Subject: [PATCH 7/7] add some comments --- examples/43-c6_coefficient.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/43-c6_coefficient.py b/examples/43-c6_coefficient.py index 235b61892..e2e05f427 100644 --- a/examples/43-c6_coefficient.py +++ b/examples/43-c6_coefficient.py @@ -26,6 +26,7 @@ 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) @@ -36,4 +37,4 @@ 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) # 45.027775533694815 \ No newline at end of file +print("Calculated C6 coefficient:", c6_val) \ No newline at end of file