Add nonadiabatic coupling vector using TDDFT-ris

examples/34-tdhf-nacv.py → 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)

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
+"""

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
+"""
+

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

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