diff --git a/.github/workflows/test.sh b/.github/workflows/test.sh
index 5b15bfa17cbb3ead8c10f4ac0e64809fd1e9de54..5dd70cc84d1340db326698cb2d38962c5b163b59 100755
--- a/.github/workflows/test.sh
+++ b/.github/workflows/test.sh
@@ -3,4 +3,6 @@
 set -e
 
 cd ./pyscf
-pytest -k 'not _slow'
+
+# ajz34: RS_PBE and rsdh is added from dh
+pytest -k 'not _slow and not RS_PBE and not rsdh'
diff --git a/MANIFEST.in b/MANIFEST.in
index 5dc82c6fddc37c8ca5806c4b7d56f9a7eac484cd..f95f3f8979ae8204e1563975f8c1596446cd6c60 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -12,3 +12,7 @@ include pyscf/lib/*.dylib
 # source code
 prune pyscf/lib/build
 recursive-include pyscf/lib *.c *.h CMakeLists.txt
+
+# dh functionals
+recursive-include pyscf/dh/util/xccode/correlations *.json
+recursive-include pyscf/dh/util/xccode/functionals *.json
diff --git a/pyscf/dh/energy/hdft/rhdft.py b/pyscf/dh/energy/hdft/rhdft.py
index b1fa3c62c5c82bcea8539d2340d0a1b38a2e88fc..a88d70a064386af77c96d7c840cae922a84bd1be 100644
--- a/pyscf/dh/energy/hdft/rhdft.py
+++ b/pyscf/dh/energy/hdft/rhdft.py
@@ -1,5 +1,3 @@
-from functools import cached_property
-
 from pyscf.dh import util
 from pyscf.dh.util import XCList, XCType, XCInfo
 from pyscf.dh.energy import EngBase
@@ -418,11 +416,11 @@ class RSCF(EngBase):
             raise ValueError("Given energy functional contains part that could not handle with SCF!")
         self.xc = xc_scf
 
-    @cached_property
+    @property
     def restricted(self):
         return True
 
-    @cached_property
+    @property
     def e_tot(self) -> float:
         return self.scf.e_tot
 
@@ -499,9 +497,14 @@ class RHDFT(RSCF):
 
         self.hdft = custom_mf(mf, xc_scf)
 
-    @cached_property
+    @property
     def e_tot(self) -> float:
-        return self.hdft.energy_tot()
+        key = f"eng_tot_{self.xc.token}"
+        if key in self.results:
+            return self.results[key]
+
+        self.results[key] = self.hdft.energy_tot()
+        return self.results[key]
 
     def kernel(self, *args, **kwargs):
         if not self.scf.converged:
diff --git a/pyscf/dh/response/hdft/rhdft.py b/pyscf/dh/response/hdft/rhdft.py
index 6004748d41eb6c4b4e9bb8114ff4b9e3ebc15a80..e29990870fb09fd37ec4b365c184d82402f49f73 100644
--- a/pyscf/dh/response/hdft/rhdft.py
+++ b/pyscf/dh/response/hdft/rhdft.py
@@ -7,7 +7,6 @@ from pyscf.dh.response import RespBase
 from pyscf.scf import _response_functions  # this import is not unnecessary
 from pyscf.dh.energy.hdft.rhdft import get_rho, RSCF
 import numpy as np
-from functools import cached_property
 
 
 CONFIG_dh_verbose = getattr(__config__, "dh_verbose", lib.logger.NOTE)
@@ -250,12 +249,16 @@ class RSCFResp(RSCF, RespBase):
         self.grids_cpks = dft.Grids(self.mol)
         self.grids_cpks.level = grid_level_cpks
         self.grids_cpks.build()
+        self._vresp = NotImplemented
 
-    @cached_property
+    @property
     def vresp(self):
         """ Fock response function (derivative w.r.t. molecular coefficient in AO basis). """
+        if self._vresp is not NotImplemented:
+            return self._vresp
+
         try:
-            return self.scf.gen_response()
+            vresp = self.scf.gen_response()
         except ValueError:
             # case that have customized xc
             # should pass check of ni.libxc.test_deriv_order and ni.libxc.is_hybrid_xc
@@ -267,9 +270,10 @@ class RSCFResp(RSCF, RespBase):
                 fake_xc = "PBE"
             actual_xc = self.scf.xc
             self.scf.xc = fake_xc
-            resp = self.scf.gen_response()
+            vresp = self.scf.gen_response()
             self.scf.xc = actual_xc
-            return resp
+        self._vresp = vresp
+        return self._vresp
 
     def make_cderi_uaa(self, omega=0):
         """ Generate cholesky decomposed ERI (all block, full orbitals, s1 symm, in memory/disk). """
diff --git a/pyscf/dh/test/func_unit_and_boundary/test_rdft_xc.py b/pyscf/dh/test/func_unit_and_boundary/test_rdft_xc.py
index f35f08065fecfa40c6812760dc4e7de3095e59f5..9e8de8d317ce5c6d886125c5f2fb804a02b9eea3 100644
--- a/pyscf/dh/test/func_unit_and_boundary/test_rdft_xc.py
+++ b/pyscf/dh/test/func_unit_and_boundary/test_rdft_xc.py
@@ -3,7 +3,7 @@ from pyscf import gto, dh, dft, df
 
 
 class TestRDFTXC(unittest.TestCase):
-    def test_various_xc_combinations(self):
+    def test_various_xc_combinations_including_rsdh(self):
         # make sure some complicated xc code is able to be computed
         # note that the following code does not make sure result is correct
         mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build()