diff --git a/examples/34-tdhf-nacv.py b/examples/34-tddft-nacv.py similarity index 55% rename from examples/34-tdhf-nacv.py rename to examples/34-tddft-nacv.py index 54dcf652..2c1aeb50 100644 --- a/examples/34-tdhf-nacv.py +++ b/examples/34-tddft-nacv.py @@ -14,7 +14,7 @@ # limitations under the License. ''' -Nonadiabatic coupling vectors between ground and excited states for RHF +Nonadiabatic coupling vectors for TDRHF and TDRKS ''' # This example will gives the derivative coupling (DC), @@ -24,6 +24,7 @@ import pyscf import gpu4pyscf from gpu4pyscf.scf import hf +from gpu4pyscf.dft import rks atom = ''' O 0.0000000000 -0.0000000000 0.1174000000 @@ -39,8 +40,8 @@ td = mf.TDA().set(nstates=5) # TDHF is OK td.kernel() # [ 9.21540892 10.99036172 11.83380819 13.62301694 15.06349085] -nac = td.NAC() -nac.state=(0,1) # same as (1,0) 0 means ground state, 1 means the first excited state +nac = td.nac_method() +nac.states=(0,1) # same as (1,0) 0 means ground state, 1 means the first excited state nac.kernel() ''' --------- TDA nonadiabatic derivative coupling for state 0 and 1---------- @@ -66,6 +67,57 @@ ---------------------------------------------- ''' +print('-----------------------------------------------------') +print("Non-adiabatic coupling matrix element (NACME) between ground and first excited state") +print(nac.de) +print('-----------------------------------------------------') +print("NACME between ground and first excited state scaled by E (/E_ex)") +print(nac.de_scaled) +print('-----------------------------------------------------') +print("NACME between ground and first excited state with ETF (electron translation factor)") +# Without including the contribution of the electron translation factor (ETF), for some molecules, +# the non-adiabatic coupling matrix element (NACME) may lack translational invariance, +# which can further lead to errors in subsequent calculations such as MD simulations. +# In this case, it is necessary to use the NACME that takes the ETF into account. +print(nac.de_etf) +print('-----------------------------------------------------') +print("NACME between ground and first excited state with ETF (electron translation factor) scaled by E (/E_ex)") +print(nac.de_etf_scaled) + + +mf = rks.RKS(mol, xc='b3lyp') # -76.4203783335521 +mf.kernel() + +td = mf.TDA().set(nstates=5) # TDHF is OK +td.kernel() # [ 7.63727447 9.47865422 10.00032863 11.95971483 14.06564139] + +nac = td.nac_method() +nac.states=(1,2) # same as (1,2) 1 means the first excited state, 2 means the second excited state +nac.kernel() +""" +--------- TDA nonadiabatic derivative coupling for states 1 and 2---------- + x y z +0 O -0.1134916788 -0.0000000000 0.0000000000 +1 H 0.0639158824 0.0000000000 0.0424664269 +2 H 0.0639158824 -0.0000000000 -0.0424664269 +--------- TDA nonadiabatic derivative coupling for states 1 and 2 after E scaled (divided by E)---------- + x y z +0 O -1.6771477392 -0.0000000000 0.0000000000 +1 H 0.9445307246 0.0000000000 0.6275567747 +2 H 0.9445307246 -0.0000000000 -0.6275567747 +--------- TDA nonadiabatic derivative coupling for states 1 and 2 with ETF---------- + x y z +0 O -0.1180169882 0.0000000000 0.0000000000 +1 H 0.0590085046 0.0000000000 0.0438508910 +2 H 0.0590085046 -0.0000000000 -0.0438508910 +--------- TDA nonadiabatic derivative coupling for states 1 and 2 with ETF after E scaled (divided by E)---------- + x y z +0 O -1.7440214738 0.0000000000 0.0000000000 +1 H 0.8720108929 0.0000000000 0.6480159906 +2 H 0.8720108929 -0.0000000000 -0.6480159906 +---------------------------------------------- +""" + print('-----------------------------------------------------') print("Non-adiabatic coupling matrix element (NACME) between ground and first excited state") print(nac.de) diff --git a/examples/37-tddft_ris_gradient.py b/examples/37-tddft_ris_gradient.py new file mode 100644 index 00000000..21647eb5 --- /dev/null +++ b/examples/37-tddft_ris_gradient.py @@ -0,0 +1,86 @@ +#!/usr/bin/env python +# 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. + +''' +Gradient for TDDFT-RIS +''' + +import pyscf +import numpy as np +import gpu4pyscf +from gpu4pyscf.dft import rks +from gpu4pyscf.tdscf import ris + +atom = ''' +O 0.0000000000 -0.0000000000 0.1174000000 +H -0.7570000000 -0.0000000000 -0.4696000000 +H 0.7570000000 0.0000000000 -0.4696000000 +''' + +mol = pyscf.M(atom=atom, basis='def2tzvp') + +mf = rks.RKS(mol, xc='pbe0') # -76.3773133945678 +mf.kernel() + +td = mf.TDA() +td.nstates=5 +td.kernel() # [ 7.81949919 9.71029362 10.13398432 12.10163229 13.93675959] (eV) + +g = td.nuc_grad_method() +g.state=1 +g.kernel() +""" +--------- TDA gradients for state 1 ---------- + x y z +0 O 0.0000000000 0.0000000000 -0.0949023769 +1 H 0.0627472634 -0.0000000000 0.0474538726 +2 H -0.0627472634 -0.0000000000 0.0474538726 +---------------------------------------------- +""" + +td_ris = ris.TDA(mf=mf.to_gpu(), nstates=5, spectra=False, single=False, gram_schmidt=True, Ktrunc=0.0) +td_ris.conv_tol = 1.0E-4 +td_ris.kernel() # [ 7.56352157 9.65590899 10.1236409 12.09873137 13.921372 ] (eV) + +g_ris = td_ris.nuc_grad_method() +g_ris.state=1 +g_ris.kernel() +""" +--------- TDA gradients for state 1 ---------- + x y z +0 O 0.0000000106 -0.0000000000 -0.0969465222 +1 H 0.0674194961 0.0000000000 0.0484759423 +2 H -0.0674195066 0.0000000000 0.0484759501 +---------------------------------------------- +""" + +print("defference for excitation energy between TDA and TDA-ris (in eV)") +print(td.e*27.21138602 - td_ris.energies.get()) +print() +""" +[0.25597762 0.05438464 0.01034341 0.00290092 0.01538759] +""" +print("defference for gradient between TDA and TDA-ris (in Hartree/Bohr)") +print(g.de - g_ris.de) +""" +[[-1.05589088e-08 -1.97977506e-15 2.04414538e-03] + [-4.67223270e-03 9.83637092e-16 -1.02206969e-03] + [ 4.67224325e-03 1.00292678e-15 -1.02207751e-03]] +""" +print("norm of the diff") +print(np.linalg.norm(g.de - g_ris.de)) +""" +0.007065933384199997 +""" \ No newline at end of file diff --git a/examples/38-tddft_ris_nacv.py b/examples/38-tddft_ris_nacv.py new file mode 100644 index 00000000..10e114fa --- /dev/null +++ b/examples/38-tddft_ris_nacv.py @@ -0,0 +1,130 @@ +#!/usr/bin/env python +# 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. + +''' +NACV for TDDFT-RIS +''' + +# This example will gives the derivative coupling (DC), +# also known as NACME (non-adiabatic coupling matrix element) +# between ground and excited states. + +import numpy as np +import pyscf +import gpu4pyscf +from gpu4pyscf.dft import rks +from gpu4pyscf.tdscf import ris + +atom = ''' +O 0.0000000000 -0.0000000000 0.1174000000 +H -0.7570000000 -0.0000000000 -0.4696000000 +H 0.7570000000 0.0000000000 -0.4696000000 +''' + +mol = pyscf.M(atom=atom, basis='def2tzvp') + +mf = rks.RKS(mol, xc='pbe0') # -76.3773133945678 +mf.conv_tol = 1e-10 +mf.kernel() + +td = mf.TDA().set(nstates=5) +td.kernel() # [ 7.81949919 9.71029362 10.13398432 12.10163229 13.93675959] + +nac = td.nac_method() +nac.states=(1,2) +nac.kernel() +""" +--------- TDA nonadiabatic derivative coupling for states 1 and 2---------- + x y z +0 O -0.0975602441 -0.0000000000 -0.0000000000 +1 H 0.0548213338 -0.0000000000 0.0360881697 +2 H 0.0548213338 0.0000000000 -0.0360881697 +--------- TDA nonadiabatic derivative coupling for states 1 and 2 after E scaled (divided by E)---------- + x y z +0 O -1.4040391809 -0.0000000000 -0.0000000000 +1 H 0.7889617464 -0.0000000000 0.5193632370 +2 H 0.7889617464 0.0000000000 -0.5193632370 +--------- TDA nonadiabatic derivative coupling for states 1 and 2 with ETF---------- + x y z +0 O -0.0965494688 -0.0000000000 -0.0000000000 +1 H 0.0482746550 -0.0000000000 0.0378920920 +2 H 0.0482746550 0.0000000000 -0.0378920920 +--------- TDA nonadiabatic derivative coupling for states 1 and 2 with ETF after E scaled (divided by E)---------- + x y z +0 O -1.3894925990 -0.0000000000 -0.0000000000 +1 H 0.6947451570 -0.0000000000 0.5453244023 +2 H 0.6947451570 0.0000000000 -0.5453244023 +---------------------------------------------- +""" + +td_ris = ris.TDA(mf=mf.to_gpu(), nstates=5, spectra=False, single=False, gram_schmidt=True, Ktrunc=0.0) +td_ris.conv_tol = 1.0E-6 +td_ris.kernel() # [ 7.56352157 9.65590898 10.12364072 12.09873136 13.921372 ] +print(td_ris.energies.get()) + +nac_ris = td_ris.nac_method() +nac_ris.states=(1,2) +nac_ris.kernel() +""" +--------- TDA nonadiabatic derivative coupling for states 1 and 2---------- + x y z +0 O 0.1009731844 0.0000000000 -0.0000000014 +1 H -0.0575662857 -0.0000000000 -0.0380920213 +2 H -0.0575662879 0.0000000000 0.0380920227 +--------- TDA nonadiabatic derivative coupling for states 1 and 2 after E scaled (divided by E)---------- + x y z +0 O 1.3131508452 0.0000000000 -0.0000000181 +1 H -0.7486464564 -0.0000000000 -0.4953846935 +2 H -0.7486464849 0.0000000000 0.4953847117 +--------- TDA nonadiabatic derivative coupling for states 1 and 2 with ETF---------- + x y z +0 O 0.1006324741 0.0000000000 -0.0000000015 +1 H -0.0503161533 -0.0000000000 -0.0395091578 +2 H -0.0503161556 0.0000000000 0.0395091592 +--------- TDA nonadiabatic derivative coupling for states 1 and 2 with ETF after E scaled (divided by E)---------- + x y z +0 O 1.3087199253 0.0000000000 -0.0000000189 +1 H -0.6543588731 -0.0000000000 -0.5138144769 +2 H -0.6543589035 0.0000000000 0.5138144958 +---------------------------------------------- +""" + +print("defference for excitation energy between TDA and TDA-ris (in eV)") +print(td.e*27.21138602 - td_ris.energies.get()) +print() +""" +[0.25597762 0.05438464 0.01034359 0.00290093 0.01538759] +""" +print("CIS derivative coupling without ETF") +print(np.abs(nac.de_scaled) - np.abs(nac_ris.de_scaled)) +print("norm of difference", np.linalg.norm(np.abs(nac.de_scaled) - np.abs(nac_ris.de_scaled))) +print() +""" +[[ 9.08883357e-02 7.20409145e-15 -1.80700507e-08] + [ 4.03152900e-02 1.15468620e-14 2.39785435e-02] + [ 4.03152615e-02 1.77621884e-16 2.39785253e-02]] + 0.11252232092869598 +""" +print("difference for CIS derivative coupling with ETF") +print(np.abs(nac.de_etf_scaled) - np.abs(nac_ris.de_etf_scaled)) +print("norm of difference", np.linalg.norm(np.abs(nac.de_etf_scaled) - np.abs(nac_ris.de_etf_scaled))) +print() +""" +[[ 8.07726737e-02 1.38040250e-14 -1.88986440e-08] + [ 4.03862838e-02 1.32300610e-14 3.15099254e-02] + [ 4.03862535e-02 -6.54111691e-16 3.15099065e-02]] +0.10849919731015174 +""" + diff --git a/gpu4pyscf/df/grad/tdrks_ris.py b/gpu4pyscf/df/grad/tdrks_ris.py new file mode 100644 index 00000000..6c63bb39 --- /dev/null +++ b/gpu4pyscf/df/grad/tdrks_ris.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. + +from gpu4pyscf.df import int3c2e, df +from gpu4pyscf.df.grad import tdrhf as tdrhf_grad_df +from gpu4pyscf.tdscf import ris +from gpu4pyscf.grad import tdrks_ris as tdrks_ris_grad +from gpu4pyscf import __config__ + +class Gradients(tdrks_ris_grad.Gradients): + from gpu4pyscf.lib.utils import to_gpu, device + + _keys = {'with_df', 'auxbasis_response'} + def __init__(self, td): + # Whether to include the response of DF auxiliary basis when computing + # nuclear gradients of J/K matrices + tdrks_ris_grad.Gradients.__init__(self, td) + + auxbasis_response = True + get_jk = tdrhf_grad_df.get_jk + + def check_sanity(self): + assert isinstance(self.base._scf, df.df_jk._DFHF) + assert isinstance(self.base, ris.TDDFT) or isinstance(self.base, ris.TDA) + + get_veff = tdrhf_grad_df.get_veff + +Grad = Gradients \ No newline at end of file diff --git a/gpu4pyscf/df/nac/__init__.py b/gpu4pyscf/df/nac/__init__.py new file mode 100644 index 00000000..35b554f8 --- /dev/null +++ b/gpu4pyscf/df/nac/__init__.py @@ -0,0 +1,15 @@ +# 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. + +from . import tdrhf, tdrks, tdrks_ris \ No newline at end of file diff --git a/gpu4pyscf/df/nac/tdrhf.py b/gpu4pyscf/df/nac/tdrhf.py new file mode 100644 index 00000000..cf331d38 --- /dev/null +++ b/gpu4pyscf/df/nac/tdrhf.py @@ -0,0 +1,38 @@ +# 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. + +from gpu4pyscf.nac import tdrhf as tdrhf_nac +from gpu4pyscf.df.grad.tdrhf import get_jk, get_veff +from gpu4pyscf.tdscf import rhf as tdrhf +from gpu4pyscf.lib import logger +from gpu4pyscf.df import df +from gpu4pyscf.lib.cupy_helper import tag_array + + +class NAC(tdrhf_nac.NAC): + from gpu4pyscf.lib.utils import to_gpu, device + + _keys = {'with_df', 'auxbasis_response'} + def __init__(self, td): + # Whether to include the response of DF auxiliary basis when computing + # nuclear gradients of J/K matrices + tdrhf_nac.NAC.__init__(self, td) + + auxbasis_response = True + get_jk = get_jk + get_veff = get_veff + + def check_sanity(self): + assert isinstance(self.base._scf, df.df_jk._DFHF) + assert isinstance(self.base, tdrhf.TDHF) or isinstance(self.base, tdrhf.TDA) diff --git a/gpu4pyscf/df/nac/tdrks.py b/gpu4pyscf/df/nac/tdrks.py new file mode 100644 index 00000000..c513d217 --- /dev/null +++ b/gpu4pyscf/df/nac/tdrks.py @@ -0,0 +1,36 @@ +# 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. + +from gpu4pyscf.nac import tdrks as tdrks_nac +from gpu4pyscf.df.nac import tdrhf as tdrhf_nac_df +from gpu4pyscf.df import df +from gpu4pyscf.tdscf import rks as tdrks + +class NAC(tdrks_nac.NAC): + from gpu4pyscf.lib.utils import to_gpu, device + + _keys = {'with_df', 'auxbasis_response'} + def __init__(self, td): + # Whether to include the response of DF auxiliary basis when computing + # nuclear gradients of J/K matrices + tdrks_nac.NAC.__init__(self, td) + + auxbasis_response = True + get_jk = tdrhf_nac_df.get_jk + + def check_sanity(self): + assert isinstance(self.base._scf, df.df_jk._DFHF) + assert isinstance(self.base, tdrks.TDDFT) or isinstance(self.base, tdrks.TDA) + + get_veff = tdrhf_nac_df.get_veff \ No newline at end of file diff --git a/gpu4pyscf/df/nac/tdrks_ris.py b/gpu4pyscf/df/nac/tdrks_ris.py new file mode 100644 index 00000000..568b2e02 --- /dev/null +++ b/gpu4pyscf/df/nac/tdrks_ris.py @@ -0,0 +1,35 @@ +# 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. + +from gpu4pyscf.nac import tdrks_ris as tdrks_ris_nac +from gpu4pyscf.df.nac import tdrhf as tdrhf_nac_df +from gpu4pyscf.df import df +from gpu4pyscf.tdscf import ris + +class NAC(tdrks_ris_nac.NAC): + from gpu4pyscf.lib.utils import to_gpu, device + + _keys = {'with_df', 'auxbasis_response'} + def __init__(self, td): + # Whether to include the response of DF auxiliary basis when computing + # nuclear gradients of J/K matrices + tdrks_ris_nac.NAC.__init__(self, td) + + auxbasis_response = True + get_veff = tdrhf_nac_df.get_veff + get_jk = tdrhf_nac_df.get_jk + + def check_sanity(self): + assert isinstance(self.base._scf, df.df_jk._DFHF) + assert isinstance(self.base, ris.TDDFT) or isinstance(self.base, ris.TDA) diff --git a/gpu4pyscf/df/tests/test_df_tddft_ris.py b/gpu4pyscf/df/tests/test_df_tddft_ris.py new file mode 100644 index 00000000..d999103c --- /dev/null +++ b/gpu4pyscf/df/tests/test_df_tddft_ris.py @@ -0,0 +1,94 @@ +# 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 unittest +import numpy as np +import cupy as cp +import pyscf +from pyscf import lib, gto, scf, dft +from gpu4pyscf import tdscf +import gpu4pyscf + +atom = """ +O 0.0000000000 0.0000000000 0.0000000000 +H 0.0000000000 -0.7570000000 0.5870000000 +H 0.0000000000 0.7570000000 0.5870000000 +""" + +bas0 = "ccpvdz" + +def setUpModule(): + global mol + mol = pyscf.M( + atom=atom, basis=bas0, max_memory=32000, output="/dev/null", verbose=1) + + +def tearDownModule(): + global mol + mol.stdout.close() + del mol + + +def diagonalize_tda(a, nroots=5): + nocc, nvir = a.shape[:2] + nov = nocc * nvir + a = a.reshape(nov, nov) + e, xy = np.linalg.eig(np.asarray(a)) + sorted_indices = np.argsort(e) + + e_sorted = e[sorted_indices] + xy_sorted = xy[:, sorted_indices] + + e_sorted_final = e_sorted[e_sorted > 1e-3] + xy_sorted = xy_sorted[:, e_sorted > 1e-3] + return e_sorted_final[:nroots], xy_sorted[:, :nroots] + + +class KnownValues(unittest.TestCase): + def test_tda_pbe_singlet(self): + mf = dft.rks.RKS(mol, xc="pbe").density_fit().to_gpu() + mf.grids.atom_grid = (99,590) + mf.conv_tol = 1.0E-12 + mf.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True, Ktrunc=0.0) + td_ris.conv_tol = 1.0E-4 + td_ris.kernel() + + a, b = td_ris.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + ref = np.array([0.25933899, 0.33439342, 0.35638257, 0.42592451, 0.51762646]) + assert np.linalg.norm(e_diag-td_ris.energies.get()/27.21138602) < 1.0E-7 + assert np.linalg.norm(e_diag-ref) < 1.0E-7 + + def test_tdaris_pbe0_singlet(self): + mf = dft.rks.RKS(mol, xc="pbe0").density_fit().to_gpu() + mf.grids.atom_grid = (99,590) + mf.conv_tol = 1.0E-12 + mf.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True, Ktrunc=0.0) + td_ris.conv_tol = 1.0E-4 + td_ris.kernel() + + a, b = td_ris.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + ref = np.array([0.28174892, 0.35852982, 0.38054425, 0.45227567, 0.5288743]) + assert np.linalg.norm(e_diag-td_ris.energies.get()/27.21138602) < 1.0E-7 + assert np.linalg.norm(e_diag-ref) < 1.0E-7 + + +if __name__ == "__main__": + print("Full Tests for density-fitting TD-RKS-ris.") + unittest.main() diff --git a/gpu4pyscf/df/tests/test_df_tddft_ris_nac.py b/gpu4pyscf/df/tests/test_df_tddft_ris_nac.py new file mode 100644 index 00000000..1c4c155c --- /dev/null +++ b/gpu4pyscf/df/tests/test_df_tddft_ris_nac.py @@ -0,0 +1,292 @@ +# 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 unittest +import numpy as np +import cupy as cp +import pyscf +from pyscf import lib, gto, scf, dft +from gpu4pyscf import tdscf, nac +import gpu4pyscf + +atom = """ +O 0.0000000000 0.0000000000 0.0000000000 +H 0.0000000000 -0.7570000000 0.5870000000 +H 0.0000000000 0.7570000000 0.5870000000 +""" + +bas0 = "def2tzvp" + +def setUpModule(): + global mol, molpbe + mol = pyscf.M( + atom=atom, basis=bas0, max_memory=32000, output="/dev/null", verbose=1) + molpbe = pyscf.M( + atom=atom, basis="ccpvdz", max_memory=32000, output="/dev/null", verbose=1) + + +def tearDownModule(): + global mol + global molpbe + mol.stdout.close() + molpbe.stdout.close() + del mol, molpbe + + +def diagonalize_tda(a, nroots=5): + nocc, nvir = a.shape[:2] + nov = nocc * nvir + a = a.reshape(nov, nov) + e, xy = np.linalg.eig(np.asarray(a)) + sorted_indices = np.argsort(e) + + e_sorted = e[sorted_indices] + xy_sorted = xy[:, sorted_indices] + + e_sorted_final = e_sorted[e_sorted > 1e-3] + xy_sorted = xy_sorted[:, e_sorted > 1e-3] + return e_sorted_final[:nroots], xy_sorted[:, :nroots] + + +class KnownValues(unittest.TestCase): + def test_nac_pbe_tdaris_singlet_vs_ref_ge(self): + mf = dft.rks.RKS(molpbe, xc="pbe").density_fit().to_gpu() + mf.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td_ris.Ktrunc = 0.0 + nac_ris = td_ris.nac_method() + + a, b = td_ris.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + nstate = 0 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ge(nac_ris, (xI, xI*0.0), e_diag[nstate]) + + ref_e = np.array([0.25933835, 0.33439277, 0.35638221, 0.42592415, 0.51762665]) + ref_de = np.array( + [[-5.90903417e-05, 7.14375913e-17, -2.26866400e-15], + [ 2.60385843e-02, 9.51909004e-16, -2.78277910e-16], + [ 2.60385843e-02, -1.05669301e-15, -1.34098234e-15],]) + ref_de_etf = np.array( + [[-1.06809321e-01, -2.75316895e-17, -1.66754809e-15], + [ 5.34045261e-02, 8.11192182e-16, -2.83779869e-16], + [ 5.34045261e-02, -9.49822786e-16, -1.37984533e-15],]) + + # compare with previous calculation resusts + assert np.linalg.norm(e_diag - ref_e) < 1.0E-8 + assert np.linalg.norm(np.abs(ana_nac[0]) - np.abs(ref_de)) < 1.0E-5 + assert np.linalg.norm(np.abs(ana_nac[2]) - np.abs(ref_de_etf)) < 1.0E-5 + + def test_nac_pbe0_tddftris_singlet_vs_ref_ge(self): + mf = dft.rks.RKS(mol, xc="pbe0").density_fit().to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + + td_ris = tdscf.ris.TDDFT(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td_ris.conv_tol = 1.0E-4 + td_ris.Ktrunc = 0.0 + td_ris.kernel() + nac_ris = td_ris.nac_method() + nac_ris.states=(1,0) + nac_ris.kernel() + + ref_de = np.array( + [[ 7.10616608e-04, 1.67751530e-17, 2.43082576e-15], + [-2.01248594e-02, -1.27707824e-15, 2.20490394e-15], + [-2.01248611e-02, 5.34947467e-16, 1.53819879e-15],]) + ref_de_etf = np.array( + [[ 1.06135079e-01, 2.11139079e-16, 2.05232306e-15], + [-5.30675522e-02, 3.70362138e-16, 1.20929274e-15], + [-5.30675585e-02, -7.92538762e-16, 5.86483960e-16],]) + + # compare with previous calculation resusts + assert np.linalg.norm(np.abs(nac_ris.de) - np.abs(ref_de)) < 1.0E-5 + assert np.linalg.norm(np.abs(nac_ris.de_etf) - np.abs(ref_de_etf)) < 1.0E-5 + + def test_nac_camb3lyp_tdaris_singlet_vs_ref_ge(self): + mf = dft.rks.RKS(mol, xc="camb3lyp").density_fit().to_gpu() + mf.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td_ris.conv_tol = 1.0E-4 + td_ris.Ktrunc = 0.0 + td_ris.kernel() + nac_ris = td_ris.nac_method() + nac_ris.states=(1,0) + nac_ris.kernel() + + ref_de = np.array( + [[-4.89770318e-04, -4.82312005e-16, 1.17207473e-15], + [ 1.90733178e-02, -8.52029551e-16, -2.58105953e-16], + [ 1.90733197e-02, 1.54498000e-15, 2.32205946e-15],]) + ref_de_etf = np.array( + [[-1.01768244e-01, -1.60217510e-16, 4.40680529e-16], + [ 5.08839991e-02, -8.82169442e-16, 4.39921371e-17], + [ 5.08840061e-02, 1.56981811e-15, 2.64021576e-15],]) + + # compare with previous calculation resusts + assert np.linalg.norm(np.abs(nac_ris.de) - np.abs(ref_de)) < 1.0E-5 + assert np.linalg.norm(np.abs(nac_ris.de_etf) - np.abs(ref_de_etf)) < 1.0E-5 + + def test_nac_pbe_tda_singlet_vs_ref_ee(self): + mf = dft.rks.RKS(molpbe, xc="pbe").density_fit().to_gpu() + mf.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td_ris.Ktrunc = 0.0 + nac_ris = td_ris.nac_method() + + a, b = td_ris.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + # excited-excited state + nstateI = 0 + nstateJ = 1 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrks_ris.get_nacv_ee(nac_ris, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + + ref_e = np.array([0.25933835, 0.33439277, 0.35638221, 0.42592415, 0.51762665]) + ref_de = np.array( + [[ 8.34605908e-17, -1.18122143e-01, -1.38959236e-14], + [ 8.91217037e-16, 6.74132293e-02, -4.46138124e-02], + [-9.93589447e-16, 6.74132293e-02, 4.46138124e-02],]) + ref_de_etf = np.array( + [[ 8.84485527e-17, -1.23821678e-01, -1.40671554e-14], + [ 8.27342873e-16, 6.19105543e-02, -4.58616923e-02], + [-9.17992771e-16, 6.19105543e-02, 4.58616923e-02],]) + + # compare with previous calculation resusts + assert np.linalg.norm(e_diag - ref_e) < 1.0E-8 + assert np.linalg.norm(np.abs(ana_nac[0]) - np.abs(ref_de)) < 1.0E-5 + assert np.linalg.norm(np.abs(ana_nac[2]) - np.abs(ref_de_etf)) < 1.0E-5 + + def test_nac_pbe_tda_singlet_fdiff(self): + """ + Compare the analytical nacv with finite difference nacv + """ + mf = dft.rks.RKS(mol, xc="pbe").density_fit().to_gpu() + mf.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True, Ktrunc=0.0) + nac_ris = td_ris.nac_method() + a, b = td_ris.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + # ground-excited state + nstate = 0 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ge(nac_ris, (xI, xI*0.0), e_diag[nstate]) + delta = 0.001 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac_ris, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 4.0E-3 + + # excited-excited state + nstateI = 1 + nstateJ = 2 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrks_ris.get_nacv_ee(nac_ris, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta=0.005 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac_ris, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta, with_ris=True) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1.0E-5 + + def test_nac_pbe0_tda_singlet_fdiff(self): + """ + Compare the analytical nacv with finite difference nacv + """ + mf = dft.rks.RKS(mol, xc="pbe0").density_fit().to_gpu() + mf.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True, Ktrunc=0.0) + nac_ris = td_ris.nac_method() + a, b = td_ris.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + # ground-excited state + nstate = 0 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ge(nac_ris, (xI, xI*0.0), e_diag[nstate]) + delta = 0.001 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac_ris, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 4.0E-3 + + # excited-excited state + nstateI = 1 + nstateJ = 2 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrks_ris.get_nacv_ee(nac_ris, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta=0.005 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac_ris, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta, with_ris=True) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1.0E-5 + + def test_nac_pbe0_tddft_singlet_vs_ref_ee(self): + mf = dft.rks.RKS(mol, xc="pbe0").density_fit().to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + + td_ris = tdscf.ris.TDDFT(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td_ris.conv_tol = 1.0E-4 + td_ris.Ktrunc = 0.0 + td_ris.kernel() + nac_ris = td_ris.nac_method() + nac_ris.states=(1,2) + nac_ris.kernel() + + ref_de = np.array( + [[-1.28630229e-16, -1.01018827e-01, -1.39860434e-09], + [-4.70672025e-16, 5.75872007e-02, -3.81043336e-02], + [-1.16131845e-16, 5.75872029e-02, 3.81043350e-02],]) + ref_de_etf = np.array( + [[-1.35556786e-16, -1.00688286e-01, -1.46281252e-09], + [-3.65451480e-16, 5.03441246e-02, -3.95286117e-02], + [-1.79485516e-16, 5.03441269e-02, 3.95286131e-02],]) + + # compare with previous calculation resusts + assert np.linalg.norm(np.abs(nac_ris.de) - np.abs(ref_de)) < 1.0E-5 + assert np.linalg.norm(np.abs(nac_ris.de_etf) - np.abs(ref_de_etf)) < 1.0E-5 + + def test_nac_camb3lyp_tddft_singlet_vs_ref_ee(self): + mf = dft.rks.RKS(mol, xc="camb3lyp").density_fit().to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + + td_ris = tdscf.ris.TDDFT(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td_ris.conv_tol = 1.0E-4 + td_ris.Ktrunc = 0.0 + td_ris.kernel() + nac_ris = td_ris.nac_method() + nac_ris.states=(1,2) + nac_ris.kernel() + + ref_de = np.array( + [[ 1.39495980e-14, -9.05064872e-02, -2.38916839e-09], + [ 6.45116755e-16, 5.26412496e-02, -3.38915479e-02], + [ 1.62215844e-15, 5.26412531e-02, 3.38915503e-02],]) + ref_de_etf = np.array( + [[ 1.40478550e-14, -9.01595513e-02, -2.48094447e-09], + [ 6.25487158e-16, 4.50797680e-02, -3.54158761e-02], + [ 1.58165609e-15, 4.50797717e-02, 3.54158785e-02],]) + + # compare with previous calculation resusts + assert np.linalg.norm(np.abs(nac_ris.de) - np.abs(ref_de)) < 1.0E-5 + assert np.linalg.norm(np.abs(nac_ris.de_etf) - np.abs(ref_de_etf)) < 1.0E-5 + + +if __name__ == "__main__": + print("Full Tests for density-fitting TD-RKS-ris nonadiabatic coupling vectors between ground and excited state.") + unittest.main() diff --git a/gpu4pyscf/df/tests/test_df_tdrhf_nac.py b/gpu4pyscf/df/tests/test_df_tdrhf_nac.py new file mode 100644 index 00000000..2bf3c16b --- /dev/null +++ b/gpu4pyscf/df/tests/test_df_tdrhf_nac.py @@ -0,0 +1,210 @@ +# 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 unittest +import numpy as np +import cupy as cp +import pyscf +from pyscf import lib, gto, scf, dft +from gpu4pyscf import tdscf, nac +import gpu4pyscf + +atom = """ +O 0.0000000000 0.0000000000 0.0000000000 +H 0.0000000000 -0.7570000000 0.5870000000 +H 0.0000000000 0.7570000000 0.5870000000 +""" + +bas0 = "cc-pvdz" + +def setUpModule(): + global mol + mol = pyscf.M( + atom=atom, basis=bas0, max_memory=32000, output="/dev/null", verbose=1) + + +def tearDownModule(): + global mol + mol.stdout.close() + del mol + + +def diagonalize_tda(a, nroots=5): + nocc, nvir = a.shape[:2] + nov = nocc * nvir + a = a.reshape(nov, nov) + e, xy = np.linalg.eig(np.asarray(a)) + sorted_indices = np.argsort(e) + + e_sorted = e[sorted_indices] + xy_sorted = xy[:, sorted_indices] + + e_sorted_final = e_sorted[e_sorted > 1e-3] + xy_sorted = xy_sorted[:, e_sorted > 1e-3] + return e_sorted_final[:nroots], xy_sorted[:, :nroots] + + +class KnownValues(unittest.TestCase): + def test_nac_tdhf_singlet_ge_vs_direct(self): + mf = scf.RHF(mol).to_gpu() + mf.kernel() + td1 = mf.TDHF().set(nstates=5) + td1.kernel() + nac1 = td1.nac_method() + nac1.states=(1,0) + nac1.kernel() + + mf = scf.RHF(mol).density_fit().to_gpu() + mf.kernel() + td2 = mf.TDHF().set(nstates=5) + td2.kernel() + nac2 = td2.nac_method() + nac2.states=(1,0) + nac2.kernel() + assert getattr(nac2.base._scf, 'with_df', None) is not None + # Compare with direct TDDFT NACV + assert abs(np.abs(nac1.de) - np.abs(nac2.de)).max() < 1e-4 + assert abs(np.abs(nac1.de_scaled) - np.abs(nac2.de_scaled)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf) - np.abs(nac2.de_etf)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf_scaled) - np.abs(nac2.de_etf_scaled)).max() < 1e-4 + + def test_nac_tda_singlet_ge_vs_direct(self): + mf = scf.RHF(mol).to_gpu() + mf.kernel() + td1 = mf.TDA().set(nstates=5) + td1.kernel() + nac1 = td1.nac_method() + nac1.states=(2,0) + nac1.kernel() + + mf = scf.RHF(mol).density_fit().to_gpu() + mf.kernel() + td2 = mf.TDA().set(nstates=5) + td2.kernel() + nac2 = td2.nac_method() + nac2.states=(2,0) + nac2.kernel() + assert getattr(nac2.base._scf, 'with_df', None) is not None + # Compare with direct TDDFT NACV + assert abs(np.abs(nac1.de) - np.abs(nac2.de)).max() < 1e-4 + assert abs(np.abs(nac1.de_scaled) - np.abs(nac2.de_scaled)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf) - np.abs(nac2.de_etf)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf_scaled) - np.abs(nac2.de_etf_scaled)).max() < 1e-4 + + def test_nac_tda_singlet_df_fdiff(self): + """ + Compare the analytical nacv with finite difference nacv + """ + mf = scf.RHF(mol).density_fit().to_gpu() + mf.kernel() + td = mf.TDA().set(nstates=5) + nac_df = td.nac_method() + + a, b = td.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + nstate = 0 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrhf.get_nacv_ge(nac_df, (xI, xI*0.0), e_diag[0]) + delta = 0.005 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac_df, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + nstateI = 0 + nstateJ = 1 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrhf.get_nacv_ee(nac_df, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta = 0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac_df, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 5.0E-5 + + + def test_nac_tda_singlet_ee_vs_direct(self): + mf = scf.RHF(mol).to_gpu() + mf.kernel() + td = mf.TDA().set(nstates=5) + td.kernel() + nac1 = td.nac_method() + nac1.states=(1,2) + nac1.kernel() + + mf1 = scf.RHF(mol).density_fit().to_gpu() + mf1.kernel() + td2 = mf1.TDA().set(nstates=5) + td2.kernel() + nac2 = td2.nac_method() + nac2.states=(1,2) + nac2.kernel() + + assert getattr(nac2.base._scf, 'with_df', None) is not None + # Compare with direct TDDFT NACV + assert abs(np.abs(nac1.de) - np.abs(nac2.de)).max() < 1e-4 + assert abs(np.abs(nac1.de_scaled) - np.abs(nac2.de_scaled)).max() < 3e-4 + assert abs(np.abs(nac1.de_etf) - np.abs(nac2.de_etf)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf_scaled) - np.abs(nac2.de_etf_scaled)).max() < 3e-4 + + nac1.states=(1,3) + nac1.kernel() + + nac2.states=(1,3) + nac2.kernel() + + assert getattr(nac2.base._scf, 'with_df', None) is not None + # Compare with direct TDDFT NACV + assert abs(np.abs(nac1.de) - np.abs(nac2.de)).max() < 1e-4 + assert abs(np.abs(nac1.de_scaled) - np.abs(nac2.de_scaled)).max() < 3e-4 + assert abs(np.abs(nac1.de_etf) - np.abs(nac2.de_etf)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf_scaled) - np.abs(nac2.de_etf_scaled)).max() < 3e-4 + + def test_nac_tdhf_singlet_ee_vs_direct(self): + mf = scf.RHF(mol).to_gpu() + mf.kernel() + td = mf.TDHF().set(nstates=5) + td.kernel() + nac1 = td.nac_method() + nac1.states=(1,2) + nac1.kernel() + + mf2 = scf.RHF(mol).density_fit().to_gpu() + mf2.kernel() + td2 = mf2.TDHF().set(nstates=5) + td2.kernel() + nac2 = td2.nac_method() + nac2.states=(1,2) + nac2.kernel() + + assert getattr(nac2.base._scf, 'with_df', None) is not None + # Compare with direct TDDFT NACV + assert abs(np.abs(nac1.de) - np.abs(nac2.de)).max() < 1e-4 + assert abs(np.abs(nac1.de_scaled) - np.abs(nac2.de_scaled)).max() < 3e-4 + assert abs(np.abs(nac1.de_etf) - np.abs(nac2.de_etf)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf_scaled) - np.abs(nac2.de_etf_scaled)).max() < 3e-4 + + nac1.states=(1,3) + nac1.kernel() + nac2.states=(1,3) + nac2.kernel() + + assert getattr(nac2.base._scf, 'with_df', None) is not None + # Compare with direct TDDFT NACV + assert abs(np.abs(nac1.de) - np.abs(nac2.de)).max() < 1e-4 + assert abs(np.abs(nac1.de_scaled) - np.abs(nac2.de_scaled)).max() < 3e-4 + assert abs(np.abs(nac1.de_etf) - np.abs(nac2.de_etf)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf_scaled) - np.abs(nac2.de_etf_scaled)).max() < 3e-4 + + +if __name__ == "__main__": + print("Full Tests for density-fitting TD-RHF nonadiabatic coupling vectors between ground and excited states.") + unittest.main() diff --git a/gpu4pyscf/df/tests/test_df_tdrks_nac.py b/gpu4pyscf/df/tests/test_df_tdrks_nac.py new file mode 100644 index 00000000..fb421796 --- /dev/null +++ b/gpu4pyscf/df/tests/test_df_tdrks_nac.py @@ -0,0 +1,269 @@ +# 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 unittest +import numpy as np +import cupy as cp +import pyscf +from pyscf import lib, gto, scf, dft +from gpu4pyscf import tdscf, nac +import gpu4pyscf + +atom = """ +O 0.0000000000 0.0000000000 0.0000000000 +H 0.0000000000 -0.7570000000 0.5870000000 +H 0.0000000000 0.7570000000 0.5870000000 +""" + +bas0 = "cc-pvdz" + +def setUpModule(): + global mol + mol = pyscf.M( + atom=atom, basis=bas0, max_memory=32000, output="/dev/null", verbose=1) + + +def tearDownModule(): + global mol + mol.stdout.close() + del mol + + +def diagonalize_tda(a, nroots=5): + nocc, nvir = a.shape[:2] + nov = nocc * nvir + a = a.reshape(nov, nov) + e, xy = np.linalg.eig(np.asarray(a)) + sorted_indices = np.argsort(e) + + e_sorted = e[sorted_indices] + xy_sorted = xy[:, sorted_indices] + + e_sorted_final = e_sorted[e_sorted > 1e-3] + xy_sorted = xy_sorted[:, e_sorted > 1e-3] + return e_sorted_final[:nroots], xy_sorted[:, :nroots] + + +class KnownValues(unittest.TestCase): + def test_nac_pbe_tddft_singlet_ge_vs_direct(self): + mf = dft.rks.RKS(mol, xc="pbe").to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td = mf.TDDFT().set(nstates=5) + td.kernel() + nac1 = td.nac_method() + nac1.states=(1,0) + nac1.kernel() + + mf = dft.rks.RKS(mol, xc="pbe").density_fit().to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td2 = mf.TDDFT().set(nstates=5) + td2.kernel() + nac2 = td2.nac_method() + nac2.states=(1,0) + nac2.kernel() + assert getattr(nac2.base._scf, 'with_df', None) is not None + # Compare with direct TDDFT NACV + assert abs(np.abs(nac1.de) - np.abs(nac2.de)).max() < 1e-4 + assert abs(np.abs(nac1.de_scaled) - np.abs(nac2.de_scaled)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf) - np.abs(nac2.de_etf)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf_scaled) - np.abs(nac2.de_etf_scaled)).max() < 1e-4 + + def test_nac_b3lyp_tddft_singlet_ge_vs_direct(self): + mf = dft.rks.RKS(mol, xc="b3lyp").to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td = mf.TDDFT().set(nstates=5) + td.kernel() + nac1 = td.nac_method() + nac1.states=(1,0) + nac1.kernel() + + mf = dft.rks.RKS(mol, xc="b3lyp").density_fit().to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td2 = mf.TDDFT().set(nstates=5) + td2.kernel() + nac2 = td2.nac_method() + nac2.states=(1,0) + nac2.kernel() + assert getattr(nac2.base._scf, 'with_df', None) is not None + # Compare with direct TDDFT NACV + assert abs(np.abs(nac1.de) - np.abs(nac2.de)).max() < 1e-4 + assert abs(np.abs(nac1.de_scaled) - np.abs(nac2.de_scaled)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf) - np.abs(nac2.de_etf)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf_scaled) - np.abs(nac2.de_etf_scaled)).max() < 1e-4 + + def test_nac_camb3lyp_tda_singlet_ge_vs_direct(self): + mf = dft.rks.RKS(mol, xc="camb3lyp").to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td = mf.TDA().set(nstates=5) + td.kernel() + nac1 = td.nac_method() + nac1.states=(1,0) + nac1.kernel() + + mf = dft.rks.RKS(mol, xc="camb3lyp").density_fit().to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td2 = mf.TDA().set(nstates=5) + td2.kernel() + nac2 = td2.nac_method() + nac2.states=(1,0) + nac2.kernel() + assert getattr(nac2.base._scf, 'with_df', None) is not None + # Compare with direct TDDFT NACV + assert abs(np.abs(nac1.de) - np.abs(nac2.de)).max() < 1e-4 + assert abs(np.abs(nac1.de_scaled) - np.abs(nac2.de_scaled)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf) - np.abs(nac2.de_etf)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf_scaled) - np.abs(nac2.de_etf_scaled)).max() < 1e-4 + + def test_nac_pbe_tda_singlet_ee_vs_direct(self): + mf = dft.rks.RKS(mol, xc="pbe").to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td = mf.TDA().set(nstates=5) + td.kernel() + nac1 = td.nac_method() + nac1.states=(1,2) + nac1.kernel() + + mf = dft.rks.RKS(mol, xc="pbe").density_fit().to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td2 = mf.TDA().set(nstates=5) + td2.kernel() + nac2 = td2.nac_method() + nac2.states=(1,2) + nac2.kernel() + # Compare with direct TDDFT NACV + assert abs(np.abs(nac1.de) - np.abs(nac2.de)).max() < 1e-4 + assert abs(np.abs(nac1.de_scaled) - np.abs(nac2.de_scaled)).max() < 5e-4 + assert abs(np.abs(nac1.de_etf) - np.abs(nac2.de_etf)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf_scaled) - np.abs(nac2.de_etf_scaled)).max() < 5e-4 + + def test_nac_pbe_tda_singlet_df_fdiff(self): + """ + Compare the analytical nacv with finite difference nacv + """ + mf = dft.rks.RKS(mol, xc="pbe").density_fit().to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td = mf.TDA().set(nstates=5) + nac1 = td.nac_method() + a, b = td.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + # ground-excited state + nstate = 0 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ge(nac1, (xI, xI*0.0), e_diag[0]) + delta = 0.001 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac1, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + # excited-excited state + nstateI = 1 + nstateJ = 2 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta = 0.001 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + def test_nac_pbe0_tda_singlet_df_fdiff(self): + """ + Compare the analytical nacv with finite difference nacv + """ + mf = dft.rks.RKS(mol, xc="pbe0").density_fit().to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td = mf.TDA().set(nstates=5) + nac1 = td.nac_method() + a, b = td.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + # ground-excited state + nstate = 0 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ge(nac1, (xI, xI*0.0), e_diag[0]) + delta = 0.001 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac1, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + # excited-excited state + nstateI = 1 + nstateJ = 2 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta = 0.001 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + def test_nac_b3lyp_tddft_singlet_ee_vs_direct(self): + mf = dft.rks.RKS(mol, xc="b3lyp").to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td = mf.TDDFT().set(nstates=5) + td.kernel() + nac1 = td.nac_method() + nac1.states=(1,2) + nac1.kernel() + + mf = dft.rks.RKS(mol, xc="b3lyp").density_fit().to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td2 = mf.TDDFT().set(nstates=5) + td2.kernel() + nac2 = td2.nac_method() + nac2.states=(1,2) + nac2.kernel() + assert abs(np.abs(nac1.de) - np.abs(nac2.de)).max() < 1e-4 + # Compare with direct TDDFT NACV + assert abs(np.abs(nac1.de_scaled) - np.abs(nac2.de_scaled)).max() < 4e-4 + assert abs(np.abs(nac1.de_etf) - np.abs(nac2.de_etf)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf_scaled) - np.abs(nac2.de_etf_scaled)).max() < 4e-4 + + def test_nac_camb3lyp_tddft_singlet_ee_vs_direct(self): + mf = dft.rks.RKS(mol, xc="camb3lyp").to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td = mf.TDDFT().set(nstates=5) + td.kernel() + nac1 = td.nac_method() + nac1.states=(1,2) + nac1.kernel() + + mf = dft.rks.RKS(mol, xc="camb3lyp").density_fit().to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td2 = mf.TDDFT().set(nstates=5) + td2.kernel() + nac2 = td2.nac_method() + nac2.states=(1,2) + nac2.kernel() + assert abs(np.abs(nac1.de) - np.abs(nac2.de)).max() < 1e-4 + # Compare with direct TDDFT NACV + assert abs(np.abs(nac1.de_scaled) - np.abs(nac2.de_scaled)).max() < 4e-4 + assert abs(np.abs(nac1.de_etf) - np.abs(nac2.de_etf)).max() < 1e-4 + assert abs(np.abs(nac1.de_etf_scaled) - np.abs(nac2.de_etf_scaled)).max() < 4e-4 + + +if __name__ == "__main__": + print("Full Tests for density-fitting TD-RKS nonadiabatic coupling vectors between ground and excited state.") + unittest.main() diff --git a/gpu4pyscf/df/tests/test_df_tdrks_ris_grad.py b/gpu4pyscf/df/tests/test_df_tdrks_ris_grad.py new file mode 100644 index 00000000..1277c105 --- /dev/null +++ b/gpu4pyscf/df/tests/test_df_tdrks_ris_grad.py @@ -0,0 +1,230 @@ +# 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 pyscf +import numpy as np +import unittest +import pytest +from pyscf import scf, dft, tdscf +import gpu4pyscf +from gpu4pyscf import scf as gpu_scf +import gpu4pyscf.tdscf.ris as ris + +atom = """ +O 0.0000000000 0.0000000000 0.0000000000 +H 0.0000000000 -0.7570000000 0.5870000000 +H 0.0000000000 0.7570000000 0.5870000000 +""" + +bas0 = "cc-pvdz" + +def diagonalize(a, b, nroots=5): + 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, xy = np.linalg.eig(np.asarray(h)) + sorted_indices = np.argsort(e) + + e_sorted = e[sorted_indices] + xy_sorted = xy[:, sorted_indices] + + e_sorted_final = e_sorted[e_sorted > 1e-3] + xy_sorted = xy_sorted[:, e_sorted > 1e-3] + return e_sorted_final[:nroots], xy_sorted[:, :nroots] + + +def diagonalize_tda(a, nroots=5): + nocc, nvir = a.shape[:2] + nov = nocc * nvir + a = a.reshape(nov, nov) + e, xy = np.linalg.eig(np.asarray(a)) + sorted_indices = np.argsort(e) + + e_sorted = e[sorted_indices] + xy_sorted = xy[:, sorted_indices] + + e_sorted_final = e_sorted[e_sorted > 1e-3] + xy_sorted = xy_sorted[:, e_sorted > 1e-3] + return e_sorted_final[:nroots], xy_sorted[:, :nroots] + + +def cal_analytic_gradient(mol, td, tdgrad, nocc, nvir, grad_elec, tda): + a, b = td.get_ab() + + if tda: + atmlst = range(mol.natm) + e_diag, xy_diag = diagonalize_tda(a) + x = xy_diag[:, 0].reshape(nocc, nvir)*np.sqrt(0.5) + de_td = grad_elec(tdgrad, (x, 0), theta=td.theta, J_fit=td.J_fit, K_fit=td.K_fit) + gradient_ana = de_td + tdgrad.grad_nuc(atmlst=atmlst) + else: + atmlst = range(mol.natm) + e_diag, xy_diag = diagonalize(a, b) + nsize = xy_diag.shape[0]//2 + norm_1 = np.linalg.norm(xy_diag[:nsize,0]) + norm_2 = np.linalg.norm(xy_diag[nsize:,0]) + x = xy_diag[:nsize,0]*np.sqrt(0.5/(norm_1**2-norm_2**2)) + y = xy_diag[nsize:,0]*np.sqrt(0.5/(norm_1**2-norm_2**2)) + x = x.reshape(nocc, nvir) + y = y.reshape(nocc, nvir) + + de_td = grad_elec(tdgrad, (x, y), theta=td.theta, J_fit=td.J_fit, K_fit=td.K_fit) + gradient_ana = de_td + tdgrad.grad_nuc(atmlst=atmlst) + + return gradient_ana + + +def setUpModule(): + global mol + mol = pyscf.M( + atom=atom, basis=bas0, max_memory=32000, output="/dev/null", verbose=1) + + +def tearDownModule(): + global mol + mol.stdout.close() + del mol + + +def benchmark_with_finite_diff( + mol_input, xc, delta=0.1, nstates=3, lindep=1.0e-12, tda=False, with_df=False): + + mol = mol_input.copy() + + if with_df: + mf = dft.RKS(mol, xc=xc).density_fit().to_gpu() + else: + mf = dft.RKS(mol, xc=xc).to_gpu() + mf.grids.level=9 + mf.grids.prune = None + mf.run() + mo_coeff = mf.mo_coeff + mo_occ = mf.mo_occ + nao, nmo = mo_coeff.shape + nocc = int((mo_occ>0).sum()) + nvir = nmo - nocc + if tda: + td = ris.TDA(mf=mf.to_gpu(), nstates=5, single=False) + else: + td = ris.TDDFT(mf=mf.to_gpu(), nstates=5, single=False) + td.conv_tol = 1.0E-8 + td.lindep=lindep + td.Ktrunc = 0.0 + td.single = False + td.nstates = nstates + tdgrad = td.nuc_grad_method() + gradient_ana = cal_analytic_gradient(mol, td, tdgrad, nocc, nvir, gpu4pyscf.grad.tdrks_ris.grad_elec, tda) + + coords = mol.atom_coords(unit='Ang')*1.0 + natm = coords.shape[0] + grad = np.zeros((natm, 3)) + for i in range(natm): + for j in range(3): + coords_new = coords*1.0 + coords_new[i, j] += delta + mol.set_geom_(coords_new, unit='Ang') + mol.build() + if with_df: + mf_add = dft.RKS(mol, xc=xc).density_fit().to_gpu() + else: + mf_add = dft.RKS(mol, xc=xc).to_gpu() + mf_add.grids.level=9 + mf_add.grids.prune = None + mf_add.run() + if tda: + td_add = ris.TDA(mf=mf_add.to_gpu(), nstates=5, single=False) + else: + td_add = ris.TDDFT(mf=mf_add.to_gpu(), nstates=5, single=False) + td_add.conv_tol = 1.0E-8 + td_add.single = False + td_add.Ktrunc = 0.0 + a, b = td_add.get_ab() + if tda: + e1 = diagonalize_tda(a)[0] + else: + e1 = diagonalize(a, b)[0] + e_add = e1[0] + mf_add.e_tot + + coords_new = coords*1.0 + coords_new[i, j] -= delta + mol.set_geom_(coords_new, unit='Ang') + mol.build() + if with_df: + mf_minus = dft.RKS(mol, xc=xc).density_fit().to_gpu() + else: + mf_minus = dft.RKS(mol, xc=xc).to_gpu() + mf_minus.grids.level=9 + mf_minus.grids.prune = None + mf_minus.run() + if tda: + td_minus = ris.TDA(mf=mf_minus.to_gpu(), nstates=5, single=False) + else: + td_minus = ris.TDDFT(mf=mf_minus.to_gpu(), nstates=5, single=False) + td_minus.conv_tol = 1.0E-8 + td_minus.single = False + td_minus.Ktrunc = 0.0 + a, b = td_minus.get_ab() + if tda: + e1 = diagonalize_tda(a)[0] + else: + e1 = diagonalize(a, b)[0] + e_minus = e1[0] + mf_minus.e_tot + + grad[i, j] = (e_add - e_minus)/(delta*2.0)*0.52917721092 + return gradient_ana, grad + + +def _check_grad(mol, xc, tol=1e-6, lindep=1.0e-12, disp=None, tda=False, with_df=False): + grad_gpu, grad = benchmark_with_finite_diff( + mol, xc, delta=0.005, nstates=5, lindep=lindep, tda=tda, with_df=with_df) + norm_diff = np.linalg.norm(grad_gpu - grad) + assert norm_diff < tol + return grad_gpu + + +class KnownValues(unittest.TestCase): + + def test_grad_b3lyp_tda_singlet_df(self): + mol = pyscf.M(atom=atom, basis='ccpvdz') + mf = dft.RKS(mol, xc='b3lyp').to_gpu() + mf.kernel() + + td = ris.TDDFT(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td.conv_tol = 1.0E-4 + td.Ktrunc = 0.0 + td.kernel() + g = td.nuc_grad_method() + g.kernel() + + mf = dft.RKS(mol, xc='b3lyp').density_fit().to_gpu() + mf.kernel() + + td2 = ris.TDDFT(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td2.conv_tol = 1.0E-4 + td2.Ktrunc = 0.0 + td2.kernel() + g2 = td2.nuc_grad_method() + g2.kernel() + + assert np.linalg.norm(g.de - g2.de) < 1.0E-4 + + def test_grad_b3lyp_tda_singlet_df_num(self): + _check_grad(mol, xc="b3lyp", tol=1e-4, tda=True, with_df=True) + + +if __name__ == "__main__": + print("Full Tests for TD-RKS RIS Gradient") + unittest.main() diff --git a/gpu4pyscf/grad/tdrhf.py b/gpu4pyscf/grad/tdrhf.py index 5d148cb8..870e45a6 100644 --- a/gpu4pyscf/grad/tdrhf.py +++ b/gpu4pyscf/grad/tdrhf.py @@ -178,12 +178,12 @@ def fvind(x): # For singlet, closed shell ground state # this term contributes the ground state contribution. dvhf = td_grad.get_veff(mol, (dmz1doo + dmz1doo.T) * 0.5 + oo0 * 2) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf # this term will remove the unused-part from PP density. dvhf = td_grad.get_veff(mol, (dmz1doo + dmz1doo.T) * 0.5) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf if singlet: j_factor=1.0 @@ -193,11 +193,11 @@ def fvind(x): # For singlet, closed shell ground state k_factor=1.0 dvhf = td_grad.get_veff(mol, (dmxpy + dmxpy.T), j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals())*2 + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())*2) dvhf_all += dvhf*2 dvhf = td_grad.get_veff(mol, (dmxmy - dmxmy.T), 0.0, k_factor, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals())*2 + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())*2) dvhf_all += dvhf*2 time1 = log.timer('2e AO integral derivatives', *time1) diff --git a/gpu4pyscf/grad/tdrks.py b/gpu4pyscf/grad/tdrks.py index c831f68b..61c336f3 100644 --- a/gpu4pyscf/grad/tdrks.py +++ b/gpu4pyscf/grad/tdrks.py @@ -215,12 +215,12 @@ def fvind(x): # this term contributes the ground state contribution. dvhf = td_grad.get_veff(mol, (dmz1doo + dmz1doo.T) * 0.5 + oo0 * 2, j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf # this term will remove the unused-part from PP density. dvhf = td_grad.get_veff(mol, (dmz1doo + dmz1doo.T) * 0.5, j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf if singlet: j_factor=1.0 @@ -228,11 +228,11 @@ def fvind(x): j_factor=0.0 dvhf = td_grad.get_veff(mol, dmxpy + dmxpy.T, j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 dvhf = td_grad.get_veff(mol, dmxmy - dmxmy.T, j_factor=0.0, k_factor=k_factor, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 if with_k and omega != 0: @@ -242,22 +242,22 @@ def fvind(x): dvhf = td_grad.get_veff(mol, (dmz1doo + dmz1doo.T) * 0.5 + oo0 * 2, j_factor=j_factor, k_factor=k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf dvhf = td_grad.get_veff(mol, (dmz1doo + dmz1doo.T) * 0.5, j_factor=j_factor, k_factor=k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_grad.get_veff(mol, dmxpy + dmxpy.T, j_factor=j_factor, k_factor=k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 dvhf = td_grad.get_veff(mol, dmxmy - dmxmy.T, j_factor=j_factor, k_factor=k_factor, omega=omega, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 time1 = log.timer('2e AO integral derivatives', *time1) fxcz1 = _contract_xc_kernel(td_grad, mf.xc, z1ao, None, False, False, True)[0] diff --git a/gpu4pyscf/grad/tdrks_ris.py b/gpu4pyscf/grad/tdrks_ris.py index 116dfff9..03c4cac4 100644 --- a/gpu4pyscf/grad/tdrks_ris.py +++ b/gpu4pyscf/grad/tdrks_ris.py @@ -40,7 +40,8 @@ def grad_elec(td_grad, x_y, theta=None, J_fit=None, K_fit=None, singlet=True, at TDDFT X and Y amplitudes. If Y is set to 0, this function computes TDA energy gradients. """ - + if td_grad.base.Ktrunc != 0.0: + raise NotImplementedError('Ktrunc or frozen method is not supported yet') if J_fit is None: J_fit = td_grad.base.J_fit if K_fit is None: @@ -69,8 +70,6 @@ def grad_elec(td_grad, x_y, theta=None, J_fit=None, K_fit=None, singlet=True, at orbo = mo_coeff[:, :nocc] if getattr(mf, 'with_solvent', None) is not None: raise NotImplementedError('With solvent is not supported yet') - if getattr(mf, 'with_df', None) is not None: - raise NotImplementedError('With density fitting is not supported yet') dvv = contract("ai,bi->ab", xpy, xpy) + contract("ai,bi->ab", xmy, xmy) # 2 T_{ab} doo = -contract("ai,aj->ij", xpy, xpy) - contract("ai,aj->ij", xmy, xmy) # 2 T_{ij} @@ -236,12 +235,12 @@ def fvind(x): # this term contributes the ground state contribution. dvhf = td_grad.get_veff(mol, (dmz1doo + dmz1doo.T) * 0.5 + oo0 * 2, j_factor=j_factor, k_factor=k_factor) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf # this term will remove the unused-part from PP density. dvhf = td_grad.get_veff(mol, (dmz1doo + dmz1doo.T) * 0.5, j_factor=j_factor, k_factor=k_factor) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf if singlet: j_factor=1.0 @@ -249,11 +248,11 @@ def fvind(x): j_factor=0.0 dvhf = get_veff_ris(mf_J, mf_K, mol, dmxpy + dmxpy.T, j_factor=j_factor, k_factor=k_factor) for k, ia in enumerate(atmlst): - extra_force[k] += get_extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(get_extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 dvhf = get_veff_ris(mf_J, mf_K, mol, dmxmy - dmxmy.T, j_factor=0.0, k_factor=k_factor, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] += get_extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(get_extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 if with_k and omega != 0: @@ -263,22 +262,22 @@ def fvind(x): dvhf = td_grad.get_veff(mol, (dmz1doo + dmz1doo.T) * 0.5 + oo0 * 2, j_factor=j_factor, k_factor=k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf dvhf = td_grad.get_veff(mol, (dmz1doo + dmz1doo.T) * 0.5, j_factor=j_factor, k_factor=k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = get_veff_ris(mf_J, mf_K, mol, dmxpy + dmxpy.T, j_factor=j_factor, k_factor=k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] += get_extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(get_extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 dvhf = get_veff_ris(mf_J, mf_K, mol, dmxmy - dmxmy.T, j_factor=j_factor, k_factor=k_factor, omega=omega, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] += get_extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(get_extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 time1 = log.timer('2e AO integral derivatives', *time1) fxcz1 = tdrks._contract_xc_kernel(td_grad, mf.xc, z1ao, None, False, False, True)[0] @@ -325,6 +324,8 @@ def kernel(self, xy=None, state=None, singlet=None, atmlst=None): state : int Excited state ID. state = 1 means the first excited state. """ + if self.base.Ktrunc != 0.0: + raise NotImplementedError('Ktrunc or frozen method is not supported yet') log = self.base.log warn_message = "TDDFT-ris gradient is still in the experimental stage, \n" +\ "and its APIs are subject to change in future releases." diff --git a/gpu4pyscf/grad/tduhf.py b/gpu4pyscf/grad/tduhf.py index e6849ee3..ff250912 100644 --- a/gpu4pyscf/grad/tduhf.py +++ b/gpu4pyscf/grad/tduhf.py @@ -227,21 +227,21 @@ def fvind(x): mol, cp.stack((((dmz1dooa + dmz1dooa.T) * 0.25 + oo0a, (dmz1doob + dmz1doob.T) * 0.25 + oo0b)))) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf # this term will remove the unused-part from PP density. dvhf = td_grad.get_veff( mol, cp.stack((((dmz1dooa + dmz1dooa.T) * 0.25, (dmz1doob + dmz1doob.T) * 0.25)))) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_grad.get_veff(mol, cp.stack((((dmxpya + dmxpya.T) * 0.5, (dmxpyb + dmxpyb.T) * 0.5)))) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 dvhf = td_grad.get_veff(mol, cp.stack((((dmxmya - dmxmya.T) * 0.5, (dmxmyb - dmxmyb.T) * 0.5))), j_factor = 0.0, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 time1 = log.timer('2e AO integral derivatives', *time1) diff --git a/gpu4pyscf/grad/tduks.py b/gpu4pyscf/grad/tduks.py index 1bcea070..301a50ea 100644 --- a/gpu4pyscf/grad/tduks.py +++ b/gpu4pyscf/grad/tduks.py @@ -280,23 +280,23 @@ def fvind(x): dvhf = td_grad.get_veff(mol, cp.stack(((dmz1dooa + dmz1dooa.T) * 0.25 + oo0a, (dmz1doob + dmz1doob.T) * 0.25 + oo0b,)), j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf # this term will remove the unused-part from PP density. dvhf = td_grad.get_veff(mol, cp.stack(((dmz1dooa + dmz1dooa.T), (dmz1doob + dmz1doob.T))) * 0.25, j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_grad.get_veff(mol, cp.stack(((dmxpya + dmxpya.T), (dmxpyb + dmxpyb.T))) * 0.5, j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 dvhf = td_grad.get_veff(mol, cp.stack(((dmxmya - dmxmya.T), (dmxmyb - dmxmyb.T))) * 0.5, j_factor=0.0, k_factor=k_factor, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 if with_k and omega != 0: @@ -308,25 +308,25 @@ def fvind(x): (dmz1doob + dmz1doob.T) * 0.25 + oo0b)), j_factor=0.0, k_factor = k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf dvhf = td_grad.get_veff(mol, cp.stack(((dmz1dooa + dmz1dooa.T) * 0.25, (dmz1doob + dmz1doob.T) * 0.25)), j_factor=0.0, k_factor = k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_grad.get_veff(mol, cp.stack(((dmxpya + dmxpya.T) * 0.5, (dmxpyb + dmxpyb.T) * 0.5)), j_factor=0.0, k_factor = k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 dvhf = td_grad.get_veff(mol, cp.stack(((dmxmya - dmxmya.T) * 0.5, (dmxmyb - dmxmyb.T) * 0.5)), j_factor=0.0, k_factor = k_factor, omega=omega, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) * 2 + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals()) * 2) dvhf_all += dvhf * 2 time1 = log.timer('2e AO integral derivatives', *time1) diff --git a/gpu4pyscf/nac/__init__.py b/gpu4pyscf/nac/__init__.py index 61a02fb6..32b10b8a 100644 --- a/gpu4pyscf/nac/__init__.py +++ b/gpu4pyscf/nac/__init__.py @@ -1,2 +1,18 @@ +# 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. + from . import tdrhf -from . import tdrks \ No newline at end of file +from . import tdrks +from . import tdrks_ris +from . import finite_diff \ No newline at end of file diff --git a/gpu4pyscf/nac/finite_diff.py b/gpu4pyscf/nac/finite_diff.py new file mode 100644 index 00000000..a8c5549d --- /dev/null +++ b/gpu4pyscf/nac/finite_diff.py @@ -0,0 +1,213 @@ +# 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 cupy as cp +import numpy as np +from gpu4pyscf import scf, dft +from gpu4pyscf.lib import logger +from gpu4pyscf.tdscf import ris + + +def diagonalize_tda(a, nroots=5): + nocc, nvir = a.shape[:2] + nov = nocc * nvir + a = a.reshape(nov, nov) + e, xy = np.linalg.eig(np.asarray(a)) + sorted_indices = np.argsort(e) + + e_sorted = e[sorted_indices] + xy_sorted = xy[:, sorted_indices] + + e_sorted_final = e_sorted[e_sorted > 1e-3] + xy_sorted = xy_sorted[:, e_sorted > 1e-3] + return e_sorted_final[:nroots], xy_sorted[:, :nroots] + + +def change_sign(s12_ao, mo_coeff_b ,mo_coeff): + mo_coeff_new = mo_coeff*1.0 + s12_mo = mo_coeff_b.T @ s12_ao @ mo_coeff_new + for i in range(s12_mo.shape[-1]): + if s12_mo[i,i] < 0.0: + mo_coeff_new[:,i] *= -1 + if s12_mo[i,i]**2 < 1.0e-10: + max_norm = -1 + for j in range(s12_mo.shape[-1]): + if s12_mo[i,j]**2 > max_norm: + max_norm = s12_mo[i,j]**2 + idx = j + mo_coeff_new[:,i] = mo_coeff[:,idx] + if mo_coeff_b[:,i].T @ s12_ao @ mo_coeff_new[:,i] < 0.0: + mo_coeff_new[:,i] *= -1 + return mo_coeff_new + + +def get_new_mol(mol, coords, delta, iatm, icart): + coords_new = coords*1.0 + coords_new[iatm, icart] += delta + mol_new = mol.copy() + mol_new.set_geom_(coords_new, unit='Ang') + return mol_new + + +def get_mf(mol, mf, s, mo_coeff): + if isinstance(mf, dft.rks.RKS): + mf_new = dft.RKS(mol) + mf_new.xc = mf.xc + if len(mf.grids.atom_grid) > 0: + mf_new.grids.atom_grid = mf.grids.atom_grid + else: + mf_new.grids.level = mf.grids.level + else: + mf_new = scf.RHF(mol) + if getattr(mf, 'with_df', None) is not None: + mf_new = mf_new.density_fit() + mf_new.conv_tol = mf.conv_tol + mf_new.conv_tol_cpscf = mf.conv_tol_cpscf + mf_new.max_cycle = mf.max_cycle + mf_new.kernel() + assert mf_new.converged + mo_coeff_new = change_sign(s, mo_coeff, mf_new.mo_coeff) + mf_new.mo_coeff = mo_coeff_new + + return mf_new + + +def get_mf_td(mol, mf, s, mo_coeff, with_ris=False): + mf_new = get_mf(mol, mf, s, mo_coeff) + if with_ris: + td_new = ris.TDA(mf=mf_new, nstates=5, spectra=False, Ktrunc = 0.0, single=False, gram_schmidt=True) + else: + td_new = mf_new.TDA() + a, b = td_new.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + return mf_new, xy_diag + + +def get_nacv_ge(td_nac, x_yI, delta=0.001, with_ris=False, singlet=True, atmlst=None, verbose=logger.INFO): + mf = td_nac.base._scf + mol = mf.mol + + coords = mol.atom_coords(unit='Ang')*1.0 + natm = coords.shape[0] + nac = np.zeros((natm, 3)) + mo_coeff = cp.asarray(mf.mo_coeff) + mo_occ = cp.asarray(mf.mo_occ) + nao, nmo = mo_coeff.shape + nocc = int((mo_occ > 0).sum()) + nvir = nmo - nocc + xI, yI = x_yI + xI = cp.asarray(xI).reshape(nocc, nvir) + if not isinstance(yI, np.ndarray) and not isinstance(yI, cp.ndarray): + yI = cp.zeros_like(xI) + yI = yI.reshape(nocc, nvir) + + gamma = np.block([[np.zeros((nocc, nocc)), xI.get()], + [(xI.T*0.0).get(), np.zeros((nvir, nvir))]]) + gamma = cp.asarray(gamma)*2 + gamma_ao = mo_coeff @ gamma @ mo_coeff.T + s = mol.intor('int1e_ovlp') + s = cp.asarray(s) + + for iatm in range(natm): + for icart in range(3): + mol_add = get_new_mol(mol, coords, delta, iatm, icart) + mf_add = get_mf(mol_add, mf, s, mo_coeff) + mol_minus = get_new_mol(mol, coords, -delta, iatm, icart) + mf_minus = get_mf(mol_minus, mf, s, mo_coeff) + + mo_diff = (mf_add.mo_coeff - mf_minus.mo_coeff)/(delta*2.0)*0.52917721092 + dpq = mo_coeff.T @ s @ mo_diff + nac[iatm, icart] = (gamma*dpq).sum() + + nac2 = np.zeros((natm, 3)) + atmlst = range(mol.natm) + offsetdic = mol.offset_nr_by_atom() + s12_deriv = mol.intor('int1e_ipovlp') + s12_deriv = cp.asarray(s12_deriv) + for k, ia in enumerate(atmlst): + shl0, shl1, p0, p1 = offsetdic[ia] + s12_deriv_tmp = s12_deriv*1.0 + ds1_tmp = s12_deriv_tmp.transpose(0,2,1) + ds1_tmp[:,:,:p0] = 0 + ds1_tmp[:,:,p1:] = 0 + nac2[k] = cp.einsum('xij,ij->x', ds1_tmp, gamma_ao).get() + return nac - nac2 + + +def get_nacv_ee(td_nac, x_yI, x_yJ, nJ, delta=0.001, with_ris=False, singlet=True, atmlst=None, verbose=logger.INFO): + mf = td_nac.base._scf + mol = mf.mol + coords = mol.atom_coords(unit='Ang')*1.0 + natm = coords.shape[0] + nac = np.zeros((natm, 3)) + nac3 = np.zeros((natm, 3)) + mo_coeff = cp.asarray(mf.mo_coeff) + mo_occ = cp.asarray(mf.mo_occ) + nao, nmo = mo_coeff.shape + nocc = int((mo_occ > 0).sum()) + nvir = nmo - nocc + xI, yI = x_yI + xJ, yJ = x_yJ + xI = cp.asarray(xI).reshape(nocc, nvir) + if not isinstance(yI, np.ndarray) and not isinstance(yI, cp.ndarray): + yI = cp.zeros_like(xI) + yI = cp.asarray(yI).reshape(nocc, nvir) + xJ = cp.asarray(xJ).reshape(nocc, nvir) + if not isinstance(yJ, np.ndarray) and not isinstance(yJ, cp.ndarray): + yJ = cp.zeros_like(xJ) + yJ = cp.asarray(yJ).reshape(nocc, nvir) + gamma = np.block([[(-xJ@xI.T).get(), np.zeros((nocc, nvir))], + [np.zeros((nvir, nocc)), (xI.T@xJ).get()]]) * 2 + gamma = cp.asarray(gamma) + gamma_ao = mo_coeff @ gamma @ mo_coeff.T + s = mol.intor('int1e_ovlp') + s = cp.asarray(s) + for iatm in range(natm): + for icart in range(3): + mol_add = get_new_mol(mol, coords, delta, iatm, icart) + mf_add, xy_diag_add = get_mf_td(mol_add, mf, s, mo_coeff, with_ris) + mol_minus = get_new_mol(mol, coords, -delta, iatm, icart) + mf_minus, xy_diag_minus = get_mf_td(mol_minus, mf, s, mo_coeff, with_ris) + + sign1 = 1.0 + sign2 = 1.0 + xJ_add = cp.asarray(xy_diag_add[:, nJ]).reshape(nocc, nvir)*cp.sqrt(0.5) + xJ_minus = cp.asarray(xy_diag_minus[:, nJ]).reshape(nocc, nvir)*cp.sqrt(0.5) + if (xJ*xJ_add).sum() < 0.0: + sign1 = -1.0 + if (xJ*xJ_minus).sum() < 0.0: + sign2 = -1.0 + + mo_diff = (mf_add.mo_coeff - mf_minus.mo_coeff)/(delta*2.0)*0.52917721092 + dpq = mo_coeff.T @ s @ mo_diff + nac[iatm, icart] = (gamma*dpq).sum() + + t_diff = (xJ_add*sign1 - xJ_minus*sign2)/(delta*2.0)*0.52917721092 + nac3[iatm, icart] = (xI*t_diff).sum()*2 # for double occupancy + + nac2 = np.zeros((natm, 3)) + atmlst = range(mol.natm) + offsetdic = mol.offset_nr_by_atom() + s12_deriv = mol.intor('int1e_ipovlp') + s12_deriv = cp.asarray(s12_deriv) + for k, ia in enumerate(atmlst): + shl0, shl1, p0, p1 = offsetdic[ia] + s12_deriv_tmp = s12_deriv*1.0 + ds1_tmp = s12_deriv_tmp.transpose(0,2,1) + ds1_tmp[:,:,:p0] = 0 + ds1_tmp[:,:,p1:] = 0 + nac2[k] = cp.einsum('xij,ij->x', ds1_tmp, gamma_ao).get() + return nac - nac2 + nac3 + \ No newline at end of file diff --git a/gpu4pyscf/nac/tdrhf.py b/gpu4pyscf/nac/tdrhf.py index 55842623..ecabe943 100644 --- a/gpu4pyscf/nac/tdrhf.py +++ b/gpu4pyscf/nac/tdrhf.py @@ -70,8 +70,6 @@ def get_nacv_ge(td_nac, x_yI, EI, singlet=True, atmlst=None, verbose=logger.INFO orbo = mo_coeff[:, :nocc] if getattr(mf, 'with_solvent', None) is not None: raise NotImplementedError('With solvent is not supported yet') - if getattr(mf, 'with_df', None) is not None: - raise NotImplementedError('With density fitting is not supported yet') xI, yI = x_yI xI = cp.asarray(xI).reshape(nocc, nvir).T @@ -203,8 +201,6 @@ def get_nacv_ee(td_nac, x_yI, x_yJ, EI, EJ, singlet=True, atmlst=None, verbose=l orbo = mo_coeff[:, :nocc] if getattr(mf, 'with_solvent', None) is not None: raise NotImplementedError('With solvent is not supported yet') - if getattr(mf, 'with_df', None) is not None: - raise NotImplementedError('With density fitting is not supported yet') xI, yI = x_yI xJ, yJ = x_yJ @@ -373,45 +369,45 @@ def fvind(x): dvhf_all = 0 dvhf = td_nac.get_veff(mol, dmz1doo + oo0) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf # minus in the next TWO terms is due to only is needed, # thus minus the contribution from same DM ({D,D}, {P,P}). dvhf = td_nac.get_veff(mol, dmz1doo) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_nac.get_veff(mol, oo0) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf j_factor=1.0 k_factor=1.0 dvhf = td_nac.get_veff(mol, (dmxpyI + dmxpyI.T + dmxpyJ + dmxpyJ.T), j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf # minus in the next TWO terms is due to only is needed, # thus minus the contribution from same DM ({R_I^S,R_I^S} and {R_J^S,R_J^S}). dvhf = td_nac.get_veff(mol, (dmxpyI + dmxpyI.T), j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf # NOTE: minus dvhf = td_nac.get_veff(mol, (dmxpyJ + dmxpyJ.T), j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_nac.get_veff(mol, (dmxmyI - dmxmyI.T + dmxmyJ - dmxmyJ.T), 0.0, k_factor, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf dvhf = td_nac.get_veff(mol, (dmxmyI - dmxmyI.T), 0.0, k_factor, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_nac.get_veff(mol, (dmxmyJ - dmxmyJ.T), 0.0, k_factor, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf delec = dh_td*2 - ds @@ -484,7 +480,7 @@ def dump_flags(self, verbose=None): ) log.info("cphf_conv_tol = %g", self.cphf_conv_tol) log.info("cphf_max_cycle = %d", self.cphf_max_cycle) - log.info("chkfile = %s", self.chkfile) + # log.info("chkfile = %s", self.chkfile) log.info(f"States ID = {self.states}") log.info("\n") return self @@ -605,4 +601,3 @@ def solvent_response(self, dm): to_gpu = lib.to_gpu -tdscf.rhf.TDA.NAC = tdscf.rhf.TDHF.NAC = lib.class_as_method(NAC) diff --git a/gpu4pyscf/nac/tdrks.py b/gpu4pyscf/nac/tdrks.py index fcdb6741..bee161a6 100644 --- a/gpu4pyscf/nac/tdrks.py +++ b/gpu4pyscf/nac/tdrks.py @@ -73,8 +73,6 @@ def get_nacv_ge(td_nac, x_yI, EI, singlet=True, atmlst=None, verbose=logger.INFO orbo = mo_coeff[:, :nocc] if getattr(mf, 'with_solvent', None) is not None: raise NotImplementedError('With solvent is not supported yet') - if getattr(mf, 'with_df', None) is not None: - raise NotImplementedError('With density fitting is not supported yet') xI, yI = x_yI xI = cp.asarray(xI).reshape(nocc, nvir).T @@ -246,8 +244,6 @@ def get_nacv_ee(td_nac, x_yI, x_yJ, EI, EJ, singlet=True, atmlst=None, verbose=l orbo = mo_coeff[:, :nocc] if getattr(mf, 'with_solvent', None) is not None: raise NotImplementedError('With solvent is not supported yet') - if getattr(mf, 'with_df', None) is not None: - raise NotImplementedError('With density fitting is not supported yet') xI, yI = x_yI xJ, yJ = x_yJ @@ -481,44 +477,44 @@ def fvind(x): dvhf_all = 0 dvhf = td_nac.get_veff(mol, dmz1doo + oo0, j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf # minus in the next TWO terms is due to only is needed, # thus minus the contribution from same DM ({D,D}, {P,P}). dvhf = td_nac.get_veff(mol, dmz1doo, j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_nac.get_veff(mol, oo0, j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_nac.get_veff(mol, (dmxpyI + dmxpyI.T + dmxpyJ + dmxpyJ.T), j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf # minus in the next TWO terms is due to only is needed, # thus minus the contribution from same DM ({R_I^S,R_I^S} and {R_J^S,R_J^S}). dvhf = td_nac.get_veff(mol, (dmxpyI + dmxpyI.T), j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf # NOTE: minus dvhf = td_nac.get_veff(mol, (dmxpyJ + dmxpyJ.T), j_factor, k_factor) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_nac.get_veff(mol, (dmxmyI - dmxmyI.T + dmxmyJ - dmxmyJ.T), 0.0, k_factor, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf dvhf = td_nac.get_veff(mol, (dmxmyI - dmxmyI.T), 0.0, k_factor, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_nac.get_veff(mol, (dmxmyJ - dmxmyJ.T), 0.0, k_factor, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf if with_k and omega != 0: @@ -526,44 +522,44 @@ def fvind(x): k_factor = alpha - hyb dvhf = td_nac.get_veff(mol, dmz1doo + oo0, j_factor, k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf # minus in the next TWO terms is due to only is needed, # thus minus the contribution from same DM ({D,D}, {P,P}). dvhf = td_nac.get_veff(mol, dmz1doo, j_factor, k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_nac.get_veff(mol, oo0, j_factor, k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_nac.get_veff(mol, (dmxpyI + dmxpyI.T + dmxpyJ + dmxpyJ.T), j_factor, k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf # minus in the next TWO terms is due to only is needed, # thus minus the contribution from same DM ({R_I^S,R_I^S} and {R_J^S,R_J^S}). dvhf = td_nac.get_veff(mol, (dmxpyI + dmxpyI.T), j_factor, k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf # NOTE: minus dvhf = td_nac.get_veff(mol, (dmxpyJ + dmxpyJ.T), j_factor, k_factor, omega=omega) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_nac.get_veff(mol, (dmxmyI - dmxmyI.T + dmxmyJ - dmxmyJ.T), 0.0, k_factor, omega=omega, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] += mf_grad.extra_force(ia, locals()) + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all += dvhf dvhf = td_nac.get_veff(mol, (dmxmyI - dmxmyI.T), 0.0, k_factor, omega=omega, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf dvhf = td_nac.get_veff(mol, (dmxmyJ - dmxmyJ.T), 0.0, k_factor, omega=omega, hermi=2) for k, ia in enumerate(atmlst): - extra_force[k] -= mf_grad.extra_force(ia, locals()) + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) dvhf_all -= dvhf fxcz1 = tdrks._contract_xc_kernel(td_nac, mf.xc, z1aoS, None, False, False, True)[0] @@ -621,4 +617,4 @@ def get_nacv_ee(self, x_yI, x_yJ, EI, EJ, singlet, atmlst=None, verbose=logger.I as_scanner = NotImplemented -tdscf.rks.TDA.NAC = tdscf.rks.TDDFT.NAC = lib.class_as_method(NAC) + diff --git a/gpu4pyscf/nac/tdrks_ris.py b/gpu4pyscf/nac/tdrks_ris.py new file mode 100644 index 00000000..85cf3639 --- /dev/null +++ b/gpu4pyscf/nac/tdrks_ris.py @@ -0,0 +1,509 @@ +# 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. +""" +Nonadiabatic derivetive coupling matrix element calculation is now in experiment. +This module is under development. +""" + +from functools import reduce +import cupy as cp +import numpy as np +from pyscf import lib +from gpu4pyscf.lib import logger +from gpu4pyscf.grad import rhf as rhf_grad +from gpu4pyscf.grad import tdrks, tdrks_ris +from gpu4pyscf.df import int3c2e +from gpu4pyscf.dft import rks +from gpu4pyscf.lib.cupy_helper import contract +from gpu4pyscf.scf import cphf +from pyscf import __config__ +from gpu4pyscf import tdscf +from gpu4pyscf.nac import tdrks as tdrks_nac +from pyscf.data.nist import HARTREE2EV +from gpu4pyscf.tdscf.ris import get_auxmol + + +def get_nacv_ee(td_nac, x_yI, x_yJ, EI, EJ, singlet=True, atmlst=None, verbose=logger.INFO): + """ + Only supports for excited-excited states. + Quadratic-response-associated terms are all neglected. + + Ref: + [1] 10.1063/1.4903986 main reference + [2] 10.1021/acs.accounts.1c00312 + [3] 10.1063/1.4885817 + + Args: + td_nac: TDNAC object + x_yI: (xI, yI) + xI and yI are the eigenvectors corresponding to the excitation and de-excitation for state I + x_yJ: (xJ, yJ) + xJ and yJ are the eigenvectors corresponding to the excitation and de-excitation for state J + EI: energy of state I + EJ: energy of state J + + Keyword args: + singlet (bool): Whether calculate singlet states. + atmlst (list): List of atoms to calculate the NAC. + verbose (int): Verbosity level. + """ + if td_nac.base.Ktrunc != 0.0: + raise NotImplementedError('Ktrunc or frozen method is not supported yet') + log = logger.new_logger(td_nac, verbose) + theta = td_nac.base.theta + J_fit = td_nac.base.J_fit + K_fit = td_nac.base.K_fit + if not singlet: + raise ValueError('TDDFT ris only supports singlet state') + mol = td_nac.mol + mf = td_nac.base._scf + mo_coeff = cp.asarray(mf.mo_coeff) + mo_energy = cp.asarray(mf.mo_energy) + mo_occ = cp.asarray(mf.mo_occ) + nao, nmo = mo_coeff.shape + nocc = int((mo_occ > 0).sum()) + nvir = nmo - nocc + orbv = mo_coeff[:, nocc:] + orbo = mo_coeff[:, :nocc] + if getattr(mf, 'with_solvent', None) is not None: + raise NotImplementedError('With solvent is not supported yet') + + xI, yI = x_yI + xJ, yJ = x_yJ + + xI = cp.asarray(xI).reshape(nocc, nvir).T + if not isinstance(yI, np.ndarray) and not isinstance(yI, cp.ndarray): + yI = cp.zeros_like(xI) + yI = cp.asarray(yI).reshape(nocc, nvir).T + xJ = cp.asarray(xJ).reshape(nocc, nvir).T + if not isinstance(yJ, np.ndarray) and not isinstance(yJ, cp.ndarray): + yJ = cp.zeros_like(xJ) + yJ = cp.asarray(yJ).reshape(nocc, nvir).T + + xpyI = (xI + yI) + xmyI = (xI - yI) + dmxpyI = reduce(cp.dot, (orbv, xpyI, orbo.T)) + dmxmyI = reduce(cp.dot, (orbv, xmyI, orbo.T)) + xpyJ = (xJ + yJ) + xmyJ = (xJ - yJ) + dmxpyJ = reduce(cp.dot, (orbv, xpyJ, orbo.T)) + dmxmyJ = reduce(cp.dot, (orbv, xmyJ, orbo.T)) + + rIJoo =-contract('ai,aj->ij', xJ, xI) - contract('ai,aj->ij', yI, yJ) + rIJvv = contract('ai,bi->ab', xI, xJ) + contract('ai,bi->ab', yJ, yI) + TIJoo = (rIJoo + rIJoo.T) * 0.5 + TIJvv = (rIJvv + rIJvv.T) * 0.5 + dmzooIJ = reduce(cp.dot, (orbo, TIJoo, orbo.T)) * 2 + dmzooIJ += reduce(cp.dot, (orbv, TIJvv, orbv.T)) * 2 + + ni = mf._numint + ni.libxc.test_deriv_order(mf.xc, 3, raise_error=True) + omega, alpha, hyb = ni.rsh_and_hybrid_coeff(mf.xc, mol.spin) + f1voI, f1ooIJ, vxc1, k1aoIJ = tdrks._contract_xc_kernel(td_nac, mf.xc, dmxpyI, dmzooIJ, True, + False, singlet) + with_k = ni.libxc.is_hybrid_xc(mf.xc) + + auxmol_J = get_auxmol(mol=mol, theta=theta, fitting_basis=J_fit) + if K_fit == J_fit and (omega == 0 or omega is None): + log.info('K uese exactly same basis as J, and they share same set of Tensors') + auxmol_K = auxmol_J + else: + log.info('K uese different basis as J') + auxmol_K = get_auxmol(mol=mol, theta=theta, fitting_basis=K_fit) + mf_J = rks.RKS(mol).density_fit() + mf_J.with_df.auxmol = auxmol_J + mf_K = rks.RKS(mol).density_fit() + mf_K.with_df.auxmol = auxmol_K + + if with_k: + vj0IJ, vk0IJ = mf.get_jk(mol, dmzooIJ, hermi=0) + vj1I = mf_J.get_j(mol, (dmxpyI + dmxpyI.T), hermi=0) + vk1I = mf_K.get_k(mol, (dmxpyI + dmxpyI.T), hermi=0) + vk2I = mf_K.get_k(mol, (dmxmyI - dmxmyI.T), hermi=0) + vj1J = mf_J.get_j(mol, (dmxpyJ + dmxpyJ.T), hermi=0) + vk1J = mf_K.get_k(mol, (dmxpyJ + dmxpyJ.T), hermi=0) + vk2J = mf_K.get_k(mol, (dmxmyJ - dmxmyJ.T), hermi=0) + + + if not isinstance(vj0IJ, cp.ndarray): + vj0IJ = cp.asarray(vj0IJ) + if not isinstance(vk0IJ, cp.ndarray): + vk0IJ = cp.asarray(vk0IJ) + if not isinstance(vj1I, cp.ndarray): + vj1I = cp.asarray(vj1I) + if not isinstance(vk1I, cp.ndarray): + vk1I = cp.asarray(vk1I) + if not isinstance(vk2I, cp.ndarray): + vk2I = cp.asarray(vk2I) + if not isinstance(vj1J, cp.ndarray): + vj1J = cp.asarray(vj1J) + if not isinstance(vk1J, cp.ndarray): + vk1J = cp.asarray(vk1J) + if not isinstance(vk2J, cp.ndarray): + vk2J = cp.asarray(vk2J) + vk0IJ *= hyb + vk1I *= hyb + vk2I *= hyb + vk1J *= hyb + vk2J *= hyb + if omega != 0: + vk0IJ_omega = mf.get_k(mol, dmzooIJ, hermi=0, omega=omega) + vk1I_omega = mf_K.get_k(mol, (dmxpyI + dmxmyI.T), hermi=0, omega=omega) + vk2I_omega = mf_K.get_k(mol, (dmxmyI - dmxmyI.T), hermi=0, omega=omega) + vk1J_omega = mf_K.get_k(mol, (dmxpyJ + dmxpyJ.T), hermi=0, omega=omega) + vk2J_omega = mf_K.get_k(mol, (dmxmyJ - dmxmyJ.T), hermi=0, omega=omega) + if not isinstance(vk0IJ, cp.ndarray): + vk0IJ = cp.asarray(vk0IJ) + if not isinstance(vk1I, cp.ndarray): + vk1I = cp.asarray(vk1I) + if not isinstance(vk2I, cp.ndarray): + vk2I = cp.asarray(vk2I) + if not isinstance(vk1J, cp.ndarray): + vk1J = cp.asarray(vk1J) + if not isinstance(vk2J, cp.ndarray): + vk2J = cp.asarray(vk2J) + vk0IJ += vk0IJ_omega * (alpha - hyb) + vk1I += vk1I_omega * (alpha - hyb) + vk2I += vk2I_omega * (alpha - hyb) + vk1J += vk1J_omega * (alpha - hyb) + vk2J += vk2J_omega * (alpha - hyb) + + veff0doo = vj0IJ * 2 - vk0IJ + f1ooIJ[0] + wvo = reduce(cp.dot, (orbv.T, veff0doo, orbo)) * 2 + veffI = vj1I * 2 - vk1I + veffI *= 0.5 + veff0mopI = reduce(cp.dot, (mo_coeff.T, veffI, mo_coeff)) + wvo -= contract("ki,ai->ak", veff0mopI[:nocc, :nocc], xpyJ) * 2 + wvo += contract("ac,ai->ci", veff0mopI[nocc:, nocc:], xpyJ) * 2 + veffJ = vj1J * 2 - vk1J + veffJ *= 0.5 + veff0mopJ = reduce(cp.dot, (mo_coeff.T, veffJ, mo_coeff)) + wvo -= contract("ki,ai->ak", veff0mopJ[:nocc, :nocc], xpyI) * 2 + wvo += contract("ac,ai->ci", veff0mopJ[nocc:, nocc:], xpyI) * 2 + veffI = -vk2I + veffI *= 0.5 + veff0momI = reduce(cp.dot, (mo_coeff.T, veffI, mo_coeff)) + wvo -= contract("ki,ai->ak", veff0momI[:nocc, :nocc], xmyJ) * 2 + wvo += contract("ac,ai->ci", veff0momI[nocc:, nocc:], xmyJ) * 2 + veffJ = -vk2J + veffJ *= 0.5 + veff0momJ = reduce(cp.dot, (mo_coeff.T, veffJ, mo_coeff)) + wvo -= contract("ki,ai->ak", veff0momJ[:nocc, :nocc], xmyI) * 2 + wvo += contract("ac,ai->ci", veff0momJ[nocc:, nocc:], xmyI) * 2 + # The up parts are according to eq. (86) and (86) in Ref. [1] + else: + vj0IJ = mf.get_j(mol, dmzooIJ, hermi=1) + vj1I = mf_J.get_j(mol, (dmxpyI + dmxpyI.T), hermi=1) + vj1J = mf_J.get_j(mol, (dmxpyJ + dmxpyJ.T), hermi=1) + if not isinstance(vj0IJ, cp.ndarray): + vj0IJ = cp.asarray(vj0IJ) + if not isinstance(vj1I, cp.ndarray): + vj1I = cp.asarray(vj1I) + if not isinstance(vj1J, cp.ndarray): + vj1J = cp.asarray(vj1J) + + veff0doo = vj0IJ * 2 + f1ooIJ[0] + wvo = reduce(cp.dot, (orbv.T, veff0doo, orbo)) * 2 + veffI = vj1I + veff0mopI = reduce(cp.dot, (mo_coeff.T, veffI, mo_coeff)) + wvo -= contract("ki,ai->ak", veff0mopI[:nocc, :nocc], xpyJ) * 2 + wvo += contract("ac,ai->ci", veff0mopI[nocc:, nocc:], xpyJ) * 2 + veffJ = vj1J + veff0mopJ = reduce(cp.dot, (mo_coeff.T, veffJ, mo_coeff)) + wvo -= contract("ki,ai->ak", veff0mopJ[:nocc, :nocc], xpyI) * 2 + wvo += contract("ac,ai->ci", veff0mopJ[nocc:, nocc:], xpyI) * 2 + veff0momI = cp.zeros((nmo, nmo)) + veff0momJ = cp.zeros((nmo, nmo)) + + vresp = mf.gen_response(singlet=None, hermi=1) + + def fvind(x): + dm = reduce(cp.dot, (orbv, x.reshape(nvir, nocc) * 2, orbo.T)) # double occupency + v1ao = vresp(dm + dm.T) + return reduce(cp.dot, (orbv.T, v1ao, orbo)).ravel() + + z1 = cphf.solve( + fvind, + mo_energy, + mo_occ, + wvo/(EJ-EI), # only one spin, negative in cphf + max_cycle=td_nac.cphf_max_cycle, + tol=td_nac.cphf_conv_tol)[0] # eq.(80) in Ref. [1] + + z1ao = reduce(cp.dot, (orbv, z1, orbo.T)) + veff = vresp((z1ao + z1ao.T)) + fock_matrix = mf.get_fock() + fock_mo = reduce(cp.dot, (mo_coeff.T, fock_matrix, mo_coeff)) + TFoo = cp.dot(TIJoo, fock_mo[:nocc,:nocc]) + TFov = cp.dot(TIJoo, fock_mo[:nocc,nocc:]) + TFvo = cp.dot(TIJvv, fock_mo[nocc:,:nocc]) + TFvv = cp.dot(TIJvv, fock_mo[nocc:,nocc:]) + + # W is calculated, eqs. (75)~(78) in Ref. [1] + # in which g_{IJ} (86) in Ref. [1] is calculated + im0 = cp.zeros((nmo, nmo)) + im0[:nocc, :nocc] = reduce(cp.dot, (orbo.T, veff0doo, orbo)) # 1st term in Eq. (81) in Ref. [1] + im0[:nocc, :nocc]+= TFoo*2.0 # 2nd term in Eq. (81) in Ref. [1] + im0[:nocc, :nocc]+= contract("ak,ai->ik", veff0mopI[nocc:, :nocc], xpyJ) # 3rd term in Eq. (81) in Ref. [1] + im0[:nocc, :nocc]+= contract("ak,ai->ik", veff0momI[nocc:, :nocc], xmyJ) # 4th term in Eq. (81) in Ref. [1] + im0[:nocc, :nocc]+= contract("ak,ai->ik", veff0mopJ[nocc:, :nocc], xpyI) # 5th term in Eq. (81) in Ref. [1] + im0[:nocc, :nocc]+= contract("ak,ai->ik", veff0momJ[nocc:, :nocc], xmyI) # 6th term in Eq. (81) in Ref. [1] + im0[:nocc, :nocc]+=rIJoo.T*(EJ-EI) # only gamma^{IJ}(II) in Eq. (29) in Ref. [1] is considered. + + im0[:nocc, nocc:] = reduce(cp.dot, (orbo.T, veff0doo, orbv)) + im0[:nocc, nocc:]+= TFov*2.0 + im0[:nocc, nocc:]+= contract("ab,ai->ib", veff0mopI[nocc:, nocc:], xpyJ) + im0[:nocc, nocc:]+= contract("ab,ai->ib", veff0momI[nocc:, nocc:], xmyJ) + im0[:nocc, nocc:]+= contract("ab,ai->ib", veff0mopJ[nocc:, nocc:], xpyI) + im0[:nocc, nocc:]+= contract("ab,ai->ib", veff0momJ[nocc:, nocc:], xmyI) + + im0[nocc:, :nocc] = TFvo*2.0 + im0[nocc:, :nocc]+= contract("ij,ai->aj", veff0mopI[:nocc, :nocc], xpyJ) + im0[nocc:, :nocc]-= contract("ij,ai->aj", veff0momI[:nocc, :nocc], xmyJ) + im0[nocc:, :nocc]+= contract("ij,ai->aj", veff0mopJ[:nocc, :nocc], xpyI) + im0[nocc:, :nocc]-= contract("ij,ai->aj", veff0momJ[:nocc, :nocc], xmyI) + + im0[nocc:, nocc:] = TFvv*2.0 + im0[nocc:, nocc:]+= contract("ib,ai->ab", veff0mopI[:nocc, nocc:], xpyJ) + im0[nocc:, nocc:]-= contract("ib,ai->ab", veff0momI[:nocc, nocc:], xmyJ) + im0[nocc:, nocc:]+= contract("ib,ai->ab", veff0mopJ[:nocc, nocc:], xpyI) + im0[nocc:, nocc:]-= contract("ib,ai->ab", veff0momJ[:nocc, nocc:], xmyI) + im0[nocc:, nocc:]+=rIJvv.T*(EJ-EI) + + im0 = im0*0.5 + im0[:nocc, :nocc]+= reduce(cp.dot, (orbo.T, veff, orbo))*(EJ-EI)*0.5 + im0[:nocc, nocc:]+= reduce(cp.dot, (orbo.T, veff, orbv))*(EJ-EI)*0.5 + im0[:nocc, nocc:]+= cp.dot(fock_mo[nocc:,nocc:],z1).T*(EJ-EI)*0.25 + im0[nocc:, :nocc]+= cp.dot(z1, fock_mo[:nocc,:nocc]*(EJ-EI))*0.25 + # 0.5 * 0.5 first is in the equation, + # second 0.5 due to z1. + # The up parts are according to eqs. (75)~(78) in Ref. [1] + # * It should be noted that, the quadratic response part is omitted! + + im0 = reduce(cp.dot, (mo_coeff, im0, mo_coeff.T))*2 + + mf_grad = td_nac.base._scf.nuc_grad_method() + s1 = mf_grad.get_ovlp(mol) + z1aoS = (z1ao + z1ao.T)*0.5* (EJ - EI) + dmz1doo = z1aoS + dmzooIJ # P + oo0 = reduce(cp.dot, (orbo, orbo.T))*2 # D + + if atmlst is None: + atmlst = range(mol.natm) + + h1 = cp.asarray(mf_grad.get_hcore(mol)) # without 1/r like terms + s1 = cp.asarray(mf_grad.get_ovlp(mol)) + dh_td = contract("xij,ij->xi", h1, dmz1doo) + ds = contract("xij,ij->xi", s1, (im0 + im0.T)) + + dh1e_td = int3c2e.get_dh1e(mol, dmz1doo) # 1/r like terms + if mol.has_ecp(): + dh1e_td += rhf_grad.get_dh1e_ecp(mol, dmz1doo) # 1/r like terms + + j_factor = 1.0 + k_factor = 0.0 + if with_k: + k_factor = hyb + + extra_force = cp.zeros((len(atmlst), 3)) + dvhf_all = 0 + dvhf = td_nac.get_veff(mol, dmz1doo + oo0, j_factor=j_factor, k_factor=k_factor) + for k, ia in enumerate(atmlst): + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) + dvhf_all += dvhf + # minus in the next TWO terms is due to only is needed, + # thus minus the contribution from same DM ({D,D}, {P,P}). + dvhf = td_nac.get_veff(mol, dmz1doo, j_factor=j_factor, k_factor=k_factor) + for k, ia in enumerate(atmlst): + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) + dvhf_all -= dvhf + dvhf = td_nac.get_veff(mol, oo0, j_factor=j_factor, k_factor=k_factor) + for k, ia in enumerate(atmlst): + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) + dvhf_all -= dvhf + dvhf = tdrks_ris.get_veff_ris(mf_J, mf_K, mol, (dmxpyI + dmxpyI.T + dmxpyJ + dmxpyJ.T), j_factor=j_factor, k_factor=k_factor) + for k, ia in enumerate(atmlst): + extra_force[k] += cp.asarray(tdrks_ris.get_extra_force(ia, locals())) + dvhf_all += dvhf + # minus in the next TWO terms is due to only is needed, + # thus minus the contribution from same DM ({R_I^S,R_I^S} and {R_J^S,R_J^S}). + dvhf = tdrks_ris.get_veff_ris(mf_J, mf_K, mol, (dmxpyI + dmxpyI.T), j_factor=j_factor, k_factor=k_factor) + for k, ia in enumerate(atmlst): + extra_force[k] -= cp.asarray(tdrks_ris.get_extra_force(ia, locals())) + dvhf_all -= dvhf # NOTE: minus + dvhf = tdrks_ris.get_veff_ris(mf_J, mf_K, mol, (dmxpyJ + dmxpyJ.T), j_factor=j_factor, k_factor=k_factor) + for k, ia in enumerate(atmlst): + extra_force[k] -= cp.asarray(tdrks_ris.get_extra_force(ia, locals())) + dvhf_all -= dvhf + dvhf = tdrks_ris.get_veff_ris(mf_J, mf_K, mol, (dmxmyI - dmxmyI.T + dmxmyJ - dmxmyJ.T), j_factor=0.0, k_factor=k_factor, hermi=2) + for k, ia in enumerate(atmlst): + extra_force[k] += cp.asarray(tdrks_ris.get_extra_force(ia, locals())) + dvhf_all += dvhf + dvhf = tdrks_ris.get_veff_ris(mf_J, mf_K, mol, (dmxmyI - dmxmyI.T), j_factor=0.0, k_factor=k_factor, hermi=2) + for k, ia in enumerate(atmlst): + extra_force[k] -= cp.asarray(tdrks_ris.get_extra_force(ia, locals())) + dvhf_all -= dvhf + dvhf = tdrks_ris.get_veff_ris(mf_J, mf_K, mol, (dmxmyJ - dmxmyJ.T), j_factor=0.0, k_factor=k_factor, hermi=2) + for k, ia in enumerate(atmlst): + extra_force[k] -= cp.asarray(tdrks_ris.get_extra_force(ia, locals())) + dvhf_all -= dvhf + + if with_k and omega != 0: + j_factor = 0.0 + k_factor = alpha - hyb + dvhf = td_nac.get_veff(mol, dmz1doo + oo0, j_factor=j_factor, k_factor=k_factor, omega=omega) + for k, ia in enumerate(atmlst): + extra_force[k] += cp.asarray(mf_grad.extra_force(ia, locals())) + dvhf_all += dvhf + # minus in the next TWO terms is due to only is needed, + # thus minus the contribution from same DM ({D,D}, {P,P}). + dvhf = td_nac.get_veff(mol, dmz1doo, j_factor=j_factor, k_factor=k_factor, omega=omega) + for k, ia in enumerate(atmlst): + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) + dvhf_all -= dvhf + dvhf = td_nac.get_veff(mol, oo0, j_factor=j_factor, k_factor=k_factor, omega=omega) + for k, ia in enumerate(atmlst): + extra_force[k] -= cp.asarray(mf_grad.extra_force(ia, locals())) + dvhf_all -= dvhf + + dvhf = tdrks_ris.get_veff_ris(mf_J, mf_K, mol, (dmxpyI + dmxpyI.T + dmxpyJ + dmxpyJ.T), j_factor=j_factor, k_factor=k_factor, omega=omega) + for k, ia in enumerate(atmlst): + extra_force[k] += cp.asarray(tdrks_ris.get_extra_force(ia, locals())) + dvhf_all += dvhf + # minus in the next TWO terms is due to only is needed, + # thus minus the contribution from same DM ({R_I^S,R_I^S} and {R_J^S,R_J^S}). + dvhf = tdrks_ris.get_veff_ris(mf_J, mf_K, mol, (dmxpyI + dmxpyI.T), j_factor=j_factor, k_factor=k_factor, omega=omega) + for k, ia in enumerate(atmlst): + extra_force[k] -= cp.asarray(tdrks_ris.get_extra_force(ia, locals())) + dvhf_all -= dvhf # NOTE: minus + dvhf = tdrks_ris.get_veff_ris(mf_J, mf_K, mol, (dmxpyJ + dmxpyJ.T), j_factor=j_factor, k_factor=k_factor, omega=omega) + for k, ia in enumerate(atmlst): + extra_force[k] -= cp.asarray(tdrks_ris.get_extra_force(ia, locals())) + dvhf_all -= dvhf + dvhf = tdrks_ris.get_veff_ris(mf_J, mf_K, mol, (dmxmyI - dmxmyI.T + dmxmyJ - dmxmyJ.T), j_factor=0.0, k_factor=k_factor, omega=omega, hermi=2) + for k, ia in enumerate(atmlst): + extra_force[k] += cp.asarray(tdrks_ris.get_extra_force(ia, locals())) + dvhf_all += dvhf + dvhf = tdrks_ris.get_veff_ris(mf_J, mf_K, mol, (dmxmyI - dmxmyI.T), j_factor=0.0, k_factor=k_factor, omega=omega, hermi=2) + for k, ia in enumerate(atmlst): + extra_force[k] -= cp.asarray(tdrks_ris.get_extra_force(ia, locals())) + dvhf_all -= dvhf + dvhf = tdrks_ris.get_veff_ris(mf_J, mf_K, mol, (dmxmyJ - dmxmyJ.T), j_factor=0.0, k_factor=k_factor, omega=omega, hermi=2) + for k, ia in enumerate(atmlst): + extra_force[k] -= cp.asarray(tdrks_ris.get_extra_force(ia, locals())) + dvhf_all -= dvhf + + fxcz1 = tdrks._contract_xc_kernel(td_nac, mf.xc, z1aoS, None, False, False, True)[0] + veff1_0 = vxc1[1:] # from in Eq. (64) in Ref.[1] + # First two terms from in Eq. (64) in Ref.[1] + # Final term from in Eq. (64) in Ref.[1] + veff1_1 = f1ooIJ[1:] + fxcz1[1:] + + delec = dh_td*2 - ds + aoslices = mol.aoslice_by_atom() + delec = cp.asarray([cp.sum(delec[:, p0:p1], axis=1) for p0, p1 in aoslices[:, 2:]]) + dveff1_0 = cp.asarray( + [contract("xpq,pq->x", veff1_0[:, p0:p1], dmz1doo[p0:p1]) for p0, p1 in aoslices[:, 2:]]) + dveff1_0 += cp.asarray([ + contract("xpq,pq->x", veff1_0[:, p0:p1].transpose(0, 2, 1), dmz1doo[:, p0:p1],) + for p0, p1 in aoslices[:, 2:]]) + dveff1_1 = cp.asarray([contract("xpq,pq->x", veff1_1[:, p0:p1], oo0[p0:p1]) for p0, p1 in aoslices[:, 2:]]) + + rIJoo_ao = reduce(cp.dot, (orbo, rIJoo, orbo.T))*2 + rIJvv_ao = reduce(cp.dot, (orbv, rIJvv, orbv.T))*2 + rIJooS_ao = reduce(cp.dot, (orbo, TIJoo, orbo.T))*2 + rIJvvS_ao = reduce(cp.dot, (orbv, TIJvv, orbv.T))*2 + ds_oo = contract("xij,ji->xi", s1, rIJoo_ao * (EJ - EI)) + ds_vv = contract("xij,ji->xi", s1, rIJvv_ao * (EJ - EI)) + ds_oo_etf = contract("xij,ji->xi", s1, rIJooS_ao * (EJ - EI)) + ds_vv_etf = contract("xij,ji->xi", s1, rIJvvS_ao * (EJ - EI)) + dsxy = cp.asarray([cp.sum(ds_oo[:, p0:p1] + ds_vv[:, p0:p1], axis=1) for p0, p1 in aoslices[:, 2:]]) + dsxy_etf = cp.asarray([cp.sum(ds_oo_etf[:, p0:p1] + ds_vv_etf[:, p0:p1], axis=1) for p0, p1 in aoslices[:, 2:]]) + + de = 2.0 * dvhf_all + extra_force + dh1e_td + delec + dveff1_0 + dveff1_1 # Eq. (64) in Ref. [1] + + de_etf = de + dsxy_etf + de += dsxy + + de = de.get() + de_etf = de_etf.get() + return de, de/(EJ - EI), de_etf, de_etf/(EJ - EI) + + +class NAC(tdrks_nac.NAC): + + @lib.with_doc(get_nacv_ee.__doc__) + def get_nacv_ee(self, x_yI, x_yJ, EI, EJ, singlet, atmlst=None, verbose=logger.INFO): + return get_nacv_ee(self, x_yI, x_yJ, EI, EJ, singlet, atmlst, verbose) + + as_scanner = NotImplemented + + def kernel(self, xy_I=None, xy_J=None, E_I=None, E_J=None, singlet=None, atmlst=None): + + logger.warn(self, "This module is under development!!") + if self.base.Ktrunc != 0.0: + raise NotImplementedError('Ktrunc or frozen method is not supported yet') + if singlet is None: + singlet = self.base.singlet + if atmlst is None: + atmlst = self.atmlst + else: + self.atmlst = atmlst + + if self.verbose >= logger.WARN: + self.check_sanity() + if self.verbose >= logger.INFO: + self.dump_flags() + + if xy_I is None or xy_J is None: + states = sorted(self.states) + nstates = len(self.base.energies) + I, J = states + if I == J: + raise ValueError("I and J should be different.") + if I < 0 or J < 0: + raise ValueError("Excited states ID should be non-negetive integers.") + elif I > nstates or J > nstates: + raise ValueError(f"Excited state exceeds the number of states {nstates}.") + elif I == 0: + logger.info(self, f"NACV between ground and excited state {J}.") + if self.base.xy[1] is not None: + xy_I = (self.base.xy[0][J-1]*np.sqrt(0.5), self.base.xy[1][J-1]*np.sqrt(0.5)) + else: + xy_I = (self.base.xy[0][J-1]*np.sqrt(0.5), self.base.xy[0][J-1]*0.0) + E_I = self.base.energies[J-1]/HARTREE2EV + E_I = float(E_I) + self.de, self.de_scaled, self.de_etf, self.de_etf_scaled \ + = self.get_nacv_ge(xy_I, E_I, singlet, atmlst, verbose=self.verbose) + self._finalize() + else: + logger.info(self, f"NACV between excited state {I} and {J}.") + if self.base.xy[1] is not None: + xy_I = (self.base.xy[0][I-1]*np.sqrt(0.5), self.base.xy[1][I-1]*np.sqrt(0.5)) + else: + xy_I = (self.base.xy[0][I-1]*np.sqrt(0.5), self.base.xy[0][I-1]*0.0) + E_I = self.base.energies[I-1]/HARTREE2EV + E_I = float(E_I) + if self.base.xy[1] is not None: + xy_J = (self.base.xy[0][J-1]*np.sqrt(0.5), self.base.xy[1][J-1]*np.sqrt(0.5)) + else: + xy_J = (self.base.xy[0][J-1]*np.sqrt(0.5), self.base.xy[0][J-1]*0.0) + E_J = self.base.energies[J-1]/HARTREE2EV + E_J = float(E_J) + self.de, self.de_scaled, self.de_etf, self.de_etf_scaled \ + = self.get_nacv_ee(xy_I, xy_J, E_I, E_J, singlet, atmlst, verbose=self.verbose) + self._finalize() + return self.de, self.de_scaled, self.de_etf, self.de_etf_scaled + + diff --git a/gpu4pyscf/nac/tests/test_tdrhf_nac_ee.py b/gpu4pyscf/nac/tests/test_tdrhf_nac_ee.py index 9a5e8f4c..d1f9e067 100644 --- a/gpu4pyscf/nac/tests/test_tdrhf_nac_ee.py +++ b/gpu4pyscf/nac/tests/test_tdrhf_nac_ee.py @@ -26,8 +26,6 @@ H 0.0000000000 0.7570000000 0.5870000000 """ -# pyscf_25 = version.parse(pyscf.__version__) <= version.parse("2.5.0") - bas0 = "cc-pvdz" def setUpModule(): @@ -42,8 +40,23 @@ def tearDownModule(): del mol +def diagonalize_tda(a, nroots=5): + nocc, nvir = a.shape[:2] + nov = nocc * nvir + a = a.reshape(nov, nov) + e, xy = np.linalg.eig(np.asarray(a)) + sorted_indices = np.argsort(e) + + e_sorted = e[sorted_indices] + xy_sorted = xy[:, sorted_indices] + + e_sorted_final = e_sorted[e_sorted > 1e-3] + xy_sorted = xy_sorted[:, e_sorted > 1e-3] + return e_sorted_final[:nroots], xy_sorted[:, :nroots] + + class KnownValues(unittest.TestCase): - def test_grad_tdhf_singlet_bdf_qchem_12(self): + def test_nac_tdhf_singlet_bdf_qchem(self): """ benchmark from both Qchem and BDF $rem @@ -123,6 +136,47 @@ def test_grad_tdhf_singlet_bdf_qchem_12(self): 1 1 1 1 1 2 noresp $end + + ==== next is qchem output ==== + --------------------------------------------------- + CIS derivative coupling without ETF + Atom X Y Z + --------------------------------------------------- + 1 -0.000000 2.322814 -0.000000 + 2 0.000000 -1.261329 0.874838 + 3 0.000000 -1.261329 -0.874838 + --------------------------------------------------- + --------------------------------------------------- + CIS Force Matrix Element + Atom X Y Z + --------------------------------------------------- + 1 -0.000000 0.155023 -0.000000 + 2 0.000000 -0.077512 0.058107 + 3 0.000000 -0.077512 -0.058107 + --------------------------------------------------- + --------------------------------------------------- + CIS derivative coupling with ETF + Atom X Y Z + --------------------------------------------------- + 1 -0.000000 2.391220 -0.000000 + 2 0.000000 -1.195610 0.896297 + 3 0.000000 -1.195610 -0.896297 + --------------------------------------------------- + ==== next is BDF output ==== + Gradient contribution from Final-NAC(R)-Escaled + 1 -0.0000000000 2.3228388812 -0.0000000000 + 2 0.0000000001 -1.2613416810 0.8748886382 + 3 -0.0000000001 -1.2613416810 -0.8748886382 + ...... + Gradient contribution from Final-NAC(S) + 1 -0.0000000000 0.1550246190 -0.0000000000 + 2 0.0000000000 -0.0775123095 0.0581102860 + 3 -0.0000000000 -0.0775123095 -0.0581102860 + ...... + Gradient contribution from Final-NAC(S)-Escaled + 1 -0.0000000000 2.3912443625 -0.0000000000 + 2 0.0000000001 -1.1956221812 0.8963472690 + 3 -0.0000000001 -1.1956221812 -0.8963472690 """ mf = scf.RHF(mol).to_gpu() mf.kernel() @@ -161,7 +215,10 @@ def test_grad_tdhf_singlet_bdf_qchem_12(self): assert abs(np.abs(nac1.de_etf) - np.abs(ref_etf_bdf)).max() < 1e-4 assert abs(np.abs(nac1.de_etf_scaled) - np.abs(ref_etf_scaled_bdf)).max() < 1e-4 - def test_grad_tda_singlet_qchem(self): + def test_nac_tda_singlet_qchem(self): + """ + Comapre with qchem + """ mf = scf.RHF(mol).to_gpu() mf.kernel() td = mf.TDA().set(nstates=5) @@ -217,6 +274,35 @@ def test_grad_tda_singlet_qchem(self): assert abs(np.abs(nac1.de_etf) - np.abs(ref_etf_qchem)).max() < 1e-4 assert abs(np.abs(nac1.de_etf_scaled) - np.abs(ref_etf_scaled_qchem)).max() < 1e-4 + def test_nac_tda_singlet_fdiff(self): + """ + Compare with finite difference + """ + mf = scf.RHF(mol).to_gpu() + mf.kernel() + td = mf.TDA().set(nstates=5) + nac1 = gpu4pyscf.nac.tdrhf.NAC(td) + a, b = td.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + nstateI = 0 + nstateJ = 1 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrhf.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta = 0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 5.0E-5 + + nstateI = 1 + nstateJ = 2 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrhf.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta = 0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 5.0E-5 + if __name__ == "__main__": print("Full Tests for TD-RHF nonadiabatic coupling vectors between excited states.") diff --git a/gpu4pyscf/nac/tests/test_tdrhf_nac_ge.py b/gpu4pyscf/nac/tests/test_tdrhf_nac_ge.py index 7924e882..d276f667 100644 --- a/gpu4pyscf/nac/tests/test_tdrhf_nac_ge.py +++ b/gpu4pyscf/nac/tests/test_tdrhf_nac_ge.py @@ -26,8 +26,6 @@ H 0.0000000000 0.7570000000 0.5870000000 """ -# pyscf_25 = version.parse(pyscf.__version__) <= version.parse("2.5.0") - bas0 = "cc-pvdz" def setUpModule(): @@ -42,8 +40,23 @@ def tearDownModule(): del mol +def diagonalize_tda(a, nroots=5): + nocc, nvir = a.shape[:2] + nov = nocc * nvir + a = a.reshape(nov, nov) + e, xy = np.linalg.eig(np.asarray(a)) + sorted_indices = np.argsort(e) + + e_sorted = e[sorted_indices] + xy_sorted = xy[:, sorted_indices] + + e_sorted_final = e_sorted[e_sorted > 1e-3] + xy_sorted = xy_sorted[:, e_sorted > 1e-3] + return e_sorted_final[:nroots], xy_sorted[:, :nroots] + + class KnownValues(unittest.TestCase): - def test_grad_tda_singlet_qchem(self): + def test_nac_tda_singlet_qchem(self): """ benchmark from qchem $rem @@ -129,7 +142,34 @@ def test_grad_tda_singlet_qchem(self): assert abs(np.abs(nac1.de_etf) - np.abs(ref_etf)).max() < 1e-4 assert abs(np.abs(nac1.de_etf_scaled) - np.abs(ref_etf_scaled)).max() < 1e-4 - def test_grad_tdhf_singlet_qchem(self): + def test_nac_tda_singlet_fdiff(self): + """ + compare with finite difference + """ + mf = scf.RHF(mol).to_gpu() + mf.kernel() + td = mf.TDA().set(nstates=5) + nac1 = gpu4pyscf.nac.tdrhf.NAC(td) + + a, b = td.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + nstate = 0 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrhf.get_nacv_ge(nac1, (xI, xI*0.0), e_diag[0]) + delta = 0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac1, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + nstate = 1 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrhf.get_nacv_ge(nac1, (xI, xI*0.0), e_diag[0]) + delta = 0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac1, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + + def test_nac_tdhf_singlet_qchem(self): """ benchmark from Qchem $rem diff --git a/gpu4pyscf/nac/tests/test_tdrks_nac_ee.py b/gpu4pyscf/nac/tests/test_tdrks_nac_ee.py index d353148b..8e733ed7 100644 --- a/gpu4pyscf/nac/tests/test_tdrks_nac_ee.py +++ b/gpu4pyscf/nac/tests/test_tdrks_nac_ee.py @@ -27,9 +27,7 @@ H 0.0000000000 0.7570000000 0.5870000000 """ -# pyscf_25 = version.parse(pyscf.__version__) <= version.parse("2.5.0") - -bas0 = "cc-pvdz" +bas0 = "ccpvdz" def setUpModule(): global mol @@ -43,8 +41,23 @@ def tearDownModule(): del mol +def diagonalize_tda(a, nroots=5): + nocc, nvir = a.shape[:2] + nov = nocc * nvir + a = a.reshape(nov, nov) + e, xy = np.linalg.eig(np.asarray(a)) + sorted_indices = np.argsort(e) + + e_sorted = e[sorted_indices] + xy_sorted = xy[:, sorted_indices] + + e_sorted_final = e_sorted[e_sorted > 1e-3] + xy_sorted = xy_sorted[:, e_sorted > 1e-3] + return e_sorted_final[:nroots], xy_sorted[:, :nroots] + + class KnownValues(unittest.TestCase): - def test_grad_pbe_tda_singlet_qchem(self): + def test_nac_pbe_tda_singlet_qchem(self): """ $rem JOBTYPE sp @@ -130,7 +143,68 @@ def test_grad_pbe_tda_singlet_qchem(self): assert abs(np.abs(nac1.de_etf) - np.abs(ref_etf)).max() < 1e-4 assert abs(np.abs(nac1.de_etf_scaled) - np.abs(ref_etf_scaled)).max() < 1e-4 - def test_grad_b3lyp_tddft_singlet_qchem(self): + def test_nac_pbe_tda_singlet_fdiff(self): + """ + compare with finite difference + """ + mf = dft.rks.RKS(mol, xc="pbe").to_gpu() + mf.kernel() + td = mf.TDA().set(nstates=5) + nac1 = td.nac_method() + a, b = td.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + nstateI = 0 + nstateJ = 1 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta = 0.005 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 4e-3 + + nstateI = 1 + nstateJ = 2 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta = 0.005 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + def test_nac_pbe0_tda_singlet_fdiff(self): + """ + compare with finite difference + """ + mf = dft.rks.RKS(mol, xc="pbe0").to_gpu() + mf.kernel() + td = mf.TDA().set(nstates=5) + nac1 = td.nac_method() + a, b = td.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + nstateI = 0 + nstateJ = 1 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta = 0.005 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-3 + + nstateI = 1 + nstateJ = 2 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta = 0.005 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac1, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + def test_nac_b3lyp_tddft_singlet_qchem(self): + """ + benchmark from qchem + """ mf = dft.rks.RKS(mol, xc="b3lyp").to_gpu() mf.grids.atom_grid = (99,590) mf.kernel() @@ -169,7 +243,10 @@ def test_grad_b3lyp_tddft_singlet_qchem(self): assert abs(np.abs(nac1.de_etf) - np.abs(ref_etf)).max() < 1e-4 assert abs(np.abs(nac1.de_etf_scaled) - np.abs(ref_etf_scaled)).max() < 1e-4 - def test_grad_camb3lyp_tddft_singlet_qchem(self): + def test_nac_camb3lyp_tddft_singlet_qchem(self): + """ + benchmark from qchem + """ mf = dft.rks.RKS(mol, xc="camb3lyp").to_gpu() mf.grids.atom_grid = (99,590) mf.kernel() diff --git a/gpu4pyscf/nac/tests/test_tdrks_nac_ge.py b/gpu4pyscf/nac/tests/test_tdrks_nac_ge.py index 8fa22cbd..af2633f7 100644 --- a/gpu4pyscf/nac/tests/test_tdrks_nac_ge.py +++ b/gpu4pyscf/nac/tests/test_tdrks_nac_ge.py @@ -26,8 +26,6 @@ H 0.0000000000 0.7570000000 0.5870000000 """ -# pyscf_25 = version.parse(pyscf.__version__) <= version.parse("2.5.0") - bas0 = "cc-pvdz" def setUpModule(): @@ -42,8 +40,23 @@ def tearDownModule(): del mol +def diagonalize_tda(a, nroots=5): + nocc, nvir = a.shape[:2] + nov = nocc * nvir + a = a.reshape(nov, nov) + e, xy = np.linalg.eig(np.asarray(a)) + sorted_indices = np.argsort(e) + + e_sorted = e[sorted_indices] + xy_sorted = xy[:, sorted_indices] + + e_sorted_final = e_sorted[e_sorted > 1e-3] + xy_sorted = xy_sorted[:, e_sorted > 1e-3] + return e_sorted_final[:nroots], xy_sorted[:, :nroots] + + class KnownValues(unittest.TestCase): - def test_grad_pbe_tda_singlet_qchem(self): + def test_nac_pbe_tda_singlet_qchem(self): """ $rem JOBTYPE sp @@ -128,9 +141,62 @@ def test_grad_pbe_tda_singlet_qchem(self): assert abs(np.abs(nac1.de_etf) - np.abs(ref_etf)).max() < 1e-4 assert abs(np.abs(nac1.de_etf_scaled) - np.abs(ref_etf_scaled)).max() < 1e-4 - - def test_grad_b3lyp_tddft_singlet_qchem(self): + def test_nac_pbe_tda_singlet_fdiff(self): + """ + compare with finite difference + """ + mf = dft.rks.RKS(mol, xc="pbe").to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td = mf.TDA().set(nstates=5) + nac1 = gpu4pyscf.nac.tdrks.NAC(td) + a, b = td.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + nstate = 0 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ge(nac1, (xI, xI*0.0), e_diag[0]) + delta = 0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac1, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + nstate = 1 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ge(nac1, (xI, xI*0.0), e_diag[0]) + delta = 0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac1, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + def test_nac_pbe0_tda_singlet_fdiff(self): + """ + compare with finite difference + """ + mf = dft.rks.RKS(mol, xc="pbe0").to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td = mf.TDA().set(nstates=5) + nac1 = gpu4pyscf.nac.tdrks.NAC(td) + a, b = td.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + nstate = 0 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ge(nac1, (xI, xI*0.0), e_diag[0]) + delta = 0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac1, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + nstate = 1 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ge(nac1, (xI, xI*0.0), e_diag[0]) + delta = 0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac1, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + def test_nac_b3lyp_tddft_singlet_qchem(self): + """ + benchmark from qchem + """ mf = dft.rks.RKS(mol, xc="b3lyp").to_gpu() mf.grids.atom_grid = (99,590) mf.kernel() @@ -169,8 +235,10 @@ def test_grad_b3lyp_tddft_singlet_qchem(self): assert abs(np.abs(nac1.de_etf) - np.abs(ref_etf)).max() < 1e-4 assert abs(np.abs(nac1.de_etf_scaled) - np.abs(ref_etf_scaled)).max() < 1e-4 - def test_grad_camb3lyp_tda_singlet_qchem(self): - + def test_nac_camb3lyp_tda_singlet_qchem(self): + """ + benchmark from qchem + """ mf = dft.rks.RKS(mol, xc="camb3lyp").to_gpu() mf.grids.atom_grid = (99,590) mf.kernel() diff --git a/gpu4pyscf/nac/tests/test_tdrks_ris_nac_ee.py b/gpu4pyscf/nac/tests/test_tdrks_ris_nac_ee.py new file mode 100644 index 00000000..1feb74e0 --- /dev/null +++ b/gpu4pyscf/nac/tests/test_tdrks_ris_nac_ee.py @@ -0,0 +1,203 @@ +# 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 unittest +import numpy as np +import cupy as cp +import pyscf +from pyscf import lib, gto, scf, dft +from gpu4pyscf import tdscf, nac +import gpu4pyscf + +atom = """ +O 0.0000000000 0.0000000000 0.0000000000 +H 0.0000000000 -0.7570000000 0.5870000000 +H 0.0000000000 0.7570000000 0.5870000000 +""" + +bas0 = "def2tzvp" + +def setUpModule(): + global mol + mol = pyscf.M( + atom=atom, basis=bas0, max_memory=32000, output="/dev/null", verbose=1) + + +def tearDownModule(): + global mol + mol.stdout.close() + del mol + + +def diagonalize_tda(a, nroots=5): + nocc, nvir = a.shape[:2] + nov = nocc * nvir + a = a.reshape(nov, nov) + e, xy = np.linalg.eig(np.asarray(a)) + sorted_indices = np.argsort(e) + + e_sorted = e[sorted_indices] + xy_sorted = xy[:, sorted_indices] + + e_sorted_final = e_sorted[e_sorted > 1e-3] + xy_sorted = xy_sorted[:, e_sorted > 1e-3] + return e_sorted_final[:nroots], xy_sorted[:, :nroots] + + +class KnownValues(unittest.TestCase): + def test_nac_pbe_tda_singlet_vs_ref(self): + mf = dft.rks.RKS(mol, xc="pbe").to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + td = mf.TDA().set(nstates=5) + td.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td_ris.conv_tol = 1.0E-4 + td_ris.Ktrunc = 0.0 + td_ris.kernel() + nac_ris = td_ris.nac_method() + nac_ris.states=(1,2) + nac_ris.kernel() + + ref_de = np.array( + [[-7.46173726e-16, 9.35902790e-02, -2.89341627e-14], + [-5.56902476e-17, -5.37437170e-02, 3.50026779e-02], + [ 7.19306347e-16, -5.37437170e-02, -3.50026779e-02],]) + ref_de_etf = np.array( + [[-6.19856849e-16, 9.26041619e-02, -2.85174872e-14], + [-2.11973474e-16, -4.63020605e-02, 3.65102194e-02], + [ 8.29062637e-16, -4.63020605e-02, -3.65102194e-02],]) + + # compare with previous calculation resusts + assert np.linalg.norm(np.abs(nac_ris.de) - np.abs(ref_de)) < 1.0E-5 + assert np.linalg.norm(np.abs(nac_ris.de_etf) - np.abs(ref_de_etf)) < 1.0E-5 + + def test_nac_pbe_tda_singlet_fdiff(self): + """ + compare with finite difference + """ + mf = dft.rks.RKS(mol, xc="pbe").to_gpu() + mf.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, Ktrunc = 0.0, spectra=False, single=False, gram_schmidt=True) + nac_ris = td_ris.nac_method() + a, b = td_ris.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + nstateI = 0 + nstateJ = 1 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrks_ris.get_nacv_ee(nac_ris, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta=0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac_ris, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta, with_ris=True) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 3.0E-3 + + nstateI = 1 + nstateJ = 2 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrks_ris.get_nacv_ee(nac_ris, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta=0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac_ris, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta, with_ris=True) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1.0E-5 + + def test_nac_pbe0_tda_singlet_fdiff(self): + """ + compare with finite difference + """ + mf = dft.rks.RKS(mol, xc="pbe0").to_gpu() + mf.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, Ktrunc = 0.0, spectra=False, single=False, gram_schmidt=True) + nac_ris = td_ris.nac_method() + a, b = td_ris.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + nstateI = 0 + nstateJ = 1 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrks_ris.get_nacv_ee(nac_ris, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta=0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac_ris, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta, with_ris=True) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 4.0E-4 + + nstateI = 1 + nstateJ = 2 + xI = xy_diag[:, nstateI]*np.sqrt(0.5) + xJ = xy_diag[:, nstateJ]*np.sqrt(0.5) + ana_nac = nac.tdrks_ris.get_nacv_ee(nac_ris, (xI, xI*0.0), (xJ, xJ*0.0), e_diag[nstateI], e_diag[nstateJ]) + delta=0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ee(nac_ris, (xI, xI*0.0), (xJ, xJ*0.0), nstateJ, delta=delta, with_ris=True) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1.0E-5 + + def test_nac_pbe0_tddft_singlet_vs_ref(self): + mf = dft.rks.RKS(mol, xc="pbe0").to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + + td_ris = tdscf.ris.TDDFT(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td_ris.conv_tol = 1.0E-4 + td_ris.Ktrunc = 0.0 + td_ris.kernel() + nac_ris = td_ris.nac_method() + nac_ris.states=(1,2) + nac_ris.kernel() + + ref_de = np.array( + [[-1.51924941e-16, -1.00952613e-01, -1.39759701e-09], + [ 1.60794931e-16, 5.75528659e-02, -3.80797499e-02], + [ 2.01916854e-16, 5.75528682e-02, 3.80797513e-02],]) + ref_de_etf = np.array( + [[-1.81973724e-16, -1.00619961e-01, -1.46179369e-09], + [ 6.25879103e-17, 5.03099621e-02, -3.95035960e-02], + [ 2.51975848e-16, 5.03099644e-02, 3.95035975e-02],]) + + # compare with previous calculation resusts + assert np.linalg.norm(np.abs(nac_ris.de) - np.abs(ref_de)) < 1.0E-5 + assert np.linalg.norm(np.abs(nac_ris.de_etf) - np.abs(ref_de_etf)) < 1.0E-5 + + + def test_nac_camb3lyp_tddft_singlet_vs_ref(self): + mf = dft.rks.RKS(mol, xc="camb3lyp").to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + + td_ris = tdscf.ris.TDDFT(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td_ris.conv_tol = 1.0E-4 + td_ris.Ktrunc = 0.0 + td_ris.kernel() + nac_ris = td_ris.nac_method() + nac_ris.states=(1,2) + nac_ris.kernel() + + ref_de = np.array( + [[ 6.09174478e-16, -9.04369401e-02, -2.38678612e-09], + [-1.25420843e-15, 5.26052262e-02, -3.38655808e-02], + [ 5.16773594e-16, 5.26052298e-02, 3.38655832e-02]]) + ref_de_etf = np.array( + [[ 5.48922915e-16, -9.00876449e-02, -2.47854031e-09], + [-1.15755555e-15, 4.50438148e-02, -3.53895964e-02], + [ 4.84556048e-16, 4.50438185e-02, 3.53895989e-02],]) + + # compare with previous calculation resusts + assert np.linalg.norm(np.abs(nac_ris.de) - np.abs(ref_de)) < 1.0E-5 + assert np.linalg.norm(np.abs(nac_ris.de_etf) - np.abs(ref_de_etf)) < 1.0E-5 + + +if __name__ == "__main__": + print("Full Tests for TD-RKS-ris nonadiabatic coupling vectors between excited states") + unittest.main() \ No newline at end of file diff --git a/gpu4pyscf/nac/tests/test_tdrks_ris_nac_ge.py b/gpu4pyscf/nac/tests/test_tdrks_ris_nac_ge.py new file mode 100644 index 00000000..9627087c --- /dev/null +++ b/gpu4pyscf/nac/tests/test_tdrks_ris_nac_ge.py @@ -0,0 +1,186 @@ +# 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 unittest +import numpy as np +import cupy as cp +import pyscf +from pyscf import lib, gto, scf, dft +from gpu4pyscf import tdscf, nac +import gpu4pyscf + +atom = """ +O 0.0000000000 0.0000000000 0.0000000000 +H 0.0000000000 -0.7570000000 0.5870000000 +H 0.0000000000 0.7570000000 0.5870000000 +""" + +# pyscf_25 = version.parse(pyscf.__version__) <= version.parse("2.5.0") + +bas0 = "def2-tzvp" + +def setUpModule(): + global mol + mol = pyscf.M( + atom=atom, basis=bas0, max_memory=32000, output="/dev/null", verbose=1) + + +def tearDownModule(): + global mol + mol.stdout.close() + del mol + +def diagonalize_tda(a, nroots=5): + nocc, nvir = a.shape[:2] + nov = nocc * nvir + a = a.reshape(nov, nov) + e, xy = np.linalg.eig(np.asarray(a)) + sorted_indices = np.argsort(e) + + e_sorted = e[sorted_indices] + xy_sorted = xy[:, sorted_indices] + + e_sorted_final = e_sorted[e_sorted > 1e-3] + xy_sorted = xy_sorted[:, e_sorted > 1e-3] + return e_sorted_final[:nroots], xy_sorted[:, :nroots] + +class KnownValues(unittest.TestCase): + def test_nac_pbe_tdaris_singlet_vs_ref(self): + mf = dft.rks.RKS(mol, xc="pbe").to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td_ris.conv_tol = 1.0E-4 + td_ris.Ktrunc = 0.0 + td_ris.kernel() + nac_ris = td_ris.nac_method() + nac_ris.states=(1,0) + nac_ris.kernel() + + ref_de = np.array( + [[ 4.38482838e-03, -8.21102914e-14, -1.69475145e-11], + [-1.88174873e-02, 9.56036995e-13, 8.47200352e-12], + [-1.88174873e-02, -8.88653417e-13, 8.37251505e-12],]) + ref_de_etf = np.array( + [[ 9.89286619e-02, -8.19838431e-14, -1.69497938e-11], + [-4.94643370e-02, 9.57436679e-13, 8.47142257e-12], + [-4.94643370e-02, -8.90160772e-13, 8.37200753e-12],]) + + # compare with previous calculation resusts + assert np.linalg.norm(np.abs(nac_ris.de) - np.abs(ref_de)) < 1.0E-5 + assert np.linalg.norm(np.abs(nac_ris.de_etf) - np.abs(ref_de_etf)) < 1.0E-5 + + def test_nac_pbe_tdaris_singlet_fdiff(self): + mf = dft.rks.RKS(mol, xc="pbe").to_gpu() + mf.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + nac_ris = td_ris.nac_method() + + a, b = td_ris.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + nstate = 0 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ge(nac_ris, (xI, xI*0.0), e_diag[0]) + delta = 0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac_ris, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + nstate = 1 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ge(nac_ris, (xI, xI*0.0), e_diag[0]) + delta = 0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac_ris, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + def test_nac_pbe0_tdaris_singlet_fdiff(self): + mf = dft.rks.RKS(mol, xc="pbe0").to_gpu() + mf.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + nac_ris = td_ris.nac_method() + + a, b = td_ris.get_ab() + e_diag, xy_diag = diagonalize_tda(a) + + nstate = 0 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ge(nac_ris, (xI, xI*0.0), e_diag[0]) + delta = 0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac_ris, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + nstate = 1 + xI = xy_diag[:, nstate]*np.sqrt(0.5) + ana_nac = nac.tdrks.get_nacv_ge(nac_ris, (xI, xI*0.0), e_diag[0]) + delta = 0.0005 + fdiff_nac = nac.finite_diff.get_nacv_ge(nac_ris, (xI, xI*0.0), delta=delta) + assert np.linalg.norm(np.abs(ana_nac[1]) - np.abs(fdiff_nac)) < 1e-5 + + def test_nac_pbe0_tddftris_singlet_vs_ref(self): + mf = dft.rks.RKS(mol, xc="pbe0").to_gpu() + mf.grids.atom_grid = (99,590) + mf.kernel() + + td_ris = tdscf.ris.TDDFT(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td_ris.conv_tol = 1.0E-4 + td_ris.Ktrunc = 0.0 + td_ris.kernel() + nac_ris = td_ris.nac_method() + nac_ris.states=(1,0) + nac_ris.kernel() + + ref_de = np.array( + [[ 7.16220831e-04, 1.01353380e-12, 1.38070626e-11], + [-2.01247331e-02, 5.83533772e-12, -6.67129184e-12], + [-2.01247348e-02, -6.83851932e-12, -7.11354404e-12],]) + ref_de_etf = np.array( + [[ 1.06105659e-01, 1.01350398e-12, 1.38098156e-11], + [-5.30528420e-02, 5.83552982e-12, -6.67197889e-12], + [-5.30528484e-02, -6.83822778e-12, -7.11429737e-12],]) + + # compare with previous calculation resusts + assert np.linalg.norm(np.abs(nac_ris.de) - np.abs(ref_de)) < 1.0E-5 + assert np.linalg.norm(np.abs(nac_ris.de_etf) - np.abs(ref_de_etf)) < 1.0E-5 + + def test_nac_camb3lyp_tdaris_singlet_vs_ref(self): + mf = dft.rks.RKS(mol, xc="camb3lyp").to_gpu() + mf.kernel() + + td_ris = tdscf.ris.TDA(mf=mf, nstates=5, spectra=False, single=False, gram_schmidt=True) + td_ris.conv_tol = 1.0E-4 + td_ris.Ktrunc = 0.0 + td_ris.kernel() + nac_ris = td_ris.nac_method() + nac_ris.states=(1,0) + nac_ris.kernel() + + ref_de = np.array( + [[-4.93284817e-04, 9.49072069e-14, 1.54699499e-11], + [ 1.90729037e-02, 5.89916730e-12, -7.71312299e-12], + [ 1.90729056e-02, -6.01789629e-12, -7.71186244e-12],]) + ref_de_etf = np.array( + [[-1.01734827e-01, 9.49409210e-14, 1.54692207e-11], + [ 5.08672907e-02, 5.90047865e-12, -7.71388350e-12], + [ 5.08672977e-02, -6.01918084e-12, -7.71263288e-12],]) + + # compare with previous calculation resusts + assert np.linalg.norm(np.abs(nac_ris.de) - np.abs(ref_de)) < 1.0E-5 + assert np.linalg.norm(np.abs(nac_ris.de_etf) - np.abs(ref_de_etf)) < 1.0E-5 + +if __name__ == "__main__": + print("Full Tests for TD-RKS-ris nonadiabatic coupling vectors between ground and excited state.") + unittest.main() diff --git a/gpu4pyscf/tdscf/rhf.py b/gpu4pyscf/tdscf/rhf.py index 18b5b85c..09975896 100644 --- a/gpu4pyscf/tdscf/rhf.py +++ b/gpu4pyscf/tdscf/rhf.py @@ -414,9 +414,10 @@ def nuc_grad_method(self): from gpu4pyscf.grad import tdrhf return tdrhf.Gradients(self) - def NAC(self): + def nac_method(self): if getattr(self._scf, 'with_df', None): - raise NotImplementedError("density fitting NAC is not supported.") + from gpu4pyscf.df.nac import tdrhf + return tdrhf.NAC(self) else: from gpu4pyscf.nac import tdrhf return tdrhf.NAC(self) diff --git a/gpu4pyscf/tdscf/ris.py b/gpu4pyscf/tdscf/ris.py index 4cebede2..d72bff2b 100644 --- a/gpu4pyscf/tdscf/ris.py +++ b/gpu4pyscf/tdscf/ris.py @@ -1127,8 +1127,20 @@ def get_3T_K(self): return T_ia_K, T_ij_K, T_ab_K def nuc_grad_method(self): - from gpu4pyscf.grad import tdrks_ris - return tdrks_ris.Gradients(self) + if getattr(self._scf, 'with_df', None) is not None: + from gpu4pyscf.df.grad import tdrks_ris + return tdrks_ris.Gradients(self) + else: + from gpu4pyscf.grad import tdrks_ris + return tdrks_ris.Gradients(self) + + def nac_method(self): + if getattr(self._scf, 'with_df', None) is not None: + from gpu4pyscf.df.nac.tdrks_ris import NAC + return NAC(self) + else: + from gpu4pyscf.nac.tdrks_ris import NAC + return NAC(self) def reset(self, mol=None): if mol is not None: diff --git a/gpu4pyscf/tdscf/rks.py b/gpu4pyscf/tdscf/rks.py index 2e991f5c..1be76207 100644 --- a/gpu4pyscf/tdscf/rks.py +++ b/gpu4pyscf/tdscf/rks.py @@ -35,9 +35,10 @@ def nuc_grad_method(self): from gpu4pyscf.grad import tdrks return tdrks.Gradients(self) - def NAC(self): + def nac_method(self): if getattr(self._scf, 'with_df', None): - raise NotImplementedError("density fitting NAC is not supported.") + from gpu4pyscf.df.nac import tdrks + return tdrks.NAC(self) else: from gpu4pyscf.nac import tdrks return tdrks.NAC(self) @@ -51,9 +52,10 @@ def nuc_grad_method(self): from gpu4pyscf.grad import tdrks return tdrks.Gradients(self) - def NAC(self): + def nac_method(self): if getattr(self._scf, 'with_df', None): - raise NotImplementedError("density fitting NAC is not supported.") + from gpu4pyscf.df.nac import tdrks + return tdrks.NAC(self) else: from gpu4pyscf.nac import tdrks return tdrks.NAC(self) diff --git a/gpu4pyscf/tdscf/uhf.py b/gpu4pyscf/tdscf/uhf.py index 54fbfce1..074abe08 100644 --- a/gpu4pyscf/tdscf/uhf.py +++ b/gpu4pyscf/tdscf/uhf.py @@ -481,6 +481,9 @@ def nuc_grad_method(self): from gpu4pyscf.grad import tduhf return tduhf.Gradients(self) + def nac_method(self): + raise NotImplementedError("Nonadiabatic coupling vector for unrestricted case is not implemented.") + def _contract_multipole(tdobj, ints, hermi=True, xy=None): if xy is None: xy = tdobj.xy mo_coeff = tdobj._scf.mo_coeff diff --git a/gpu4pyscf/tdscf/uks.py b/gpu4pyscf/tdscf/uks.py index 9970630a..c808fd9d 100644 --- a/gpu4pyscf/tdscf/uks.py +++ b/gpu4pyscf/tdscf/uks.py @@ -37,6 +37,9 @@ def nuc_grad_method(self): from gpu4pyscf.grad import tduks return tduks.Gradients(self) + def nac_method(self): + raise NotImplementedError("Nonadiabatic coupling vector for unrestricted case is not implemented.") + class TDDFT(tdhf_gpu.TDHF): def nuc_grad_method(self): if getattr(self._scf, 'with_df', None): @@ -46,6 +49,9 @@ def nuc_grad_method(self): from gpu4pyscf.grad import tduks return tduks.Gradients(self) + def nac_method(self): + raise NotImplementedError("Nonadiabatic coupling vector for unrestricted case is not implemented.") + TDUKS = TDDFT SpinFlipTDA = tdhf_gpu.SpinFlipTDA SpinFlipTDDFT = tdhf_gpu.SpinFlipTDHF