diff --git a/CITATION.cff b/CITATION.cff
index e0fe04ac..29cb8162 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -8,6 +8,6 @@ authors:
 - family-names: Booth
   given-names: George H.
 title: "ebcc: Coupled cluster calculations on electron-boson systems"
-version: 1.4.5
+version: 1.5.0
 date-released: 2024-07-31
 url: "https://github.com/BoothGroup/ebcc"
diff --git a/ebcc/__init__.py b/ebcc/__init__.py
index a8576cb6..0e8a8052 100644
--- a/ebcc/__init__.py
+++ b/ebcc/__init__.py
@@ -37,7 +37,7 @@
 from __future__ import annotations
 
 """Version of the package."""
-__version__ = "1.4.5"
+__version__ = "1.5.0"
 
 """List of supported ansatz types."""
 METHOD_TYPES = ["MP", "CC", "LCC", "QCI", "QCC", "DC"]
diff --git a/ebcc/cc/base.py b/ebcc/cc/base.py
index eb4753d4..be92c5ee 100644
--- a/ebcc/cc/base.py
+++ b/ebcc/cc/base.py
@@ -782,11 +782,34 @@ def make_eb_coup_rdm(
         """
         pass
 
+    def hbar_matvec_ip_intermediates(
+        self,
+        eris: Optional[ERIsInputType] = None,
+        amplitudes: Optional[Namespace[SpinArrayType]] = None,
+    ) -> Namespace[NDArray[float]]:
+        """Compute intermediates for the IP-EOM Hamiltonian.
+
+        Args:
+            eris: Electron repulsion integrals.
+            amplitudes: Cluster amplitudes.
+
+        Returns:
+            Intermediate products for the IP-EOM Hamiltonian.
+        """
+        func, kwargs = self._load_function(
+            "hbar_matvec_ip_intermediates",
+            eris=eris,
+            amplitudes=amplitudes,
+        )
+        res: Namespace[NDArray[float]] = util.Namespace(**func(**kwargs))
+        return res
+
     def hbar_matvec_ip(
         self,
         excitations: Namespace[SpinArrayType],
         eris: Optional[ERIsInputType] = None,
         amplitudes: Optional[Namespace[SpinArrayType]] = None,
+        ints: Optional[Namespace[NDArray[float]]] = None,
     ) -> Namespace[SpinArrayType]:
         """Compute the product between a state vector and the IP-EOM Hamiltonian.
 
@@ -794,26 +817,53 @@ def hbar_matvec_ip(
             excitations: State vector as a set of excitation amplitudes.
             eris: Electron repulsion integrals.
             amplitudes: Cluster amplitudes.
+            ints: Intermediate products.
 
         Returns:
             Products between the state vectors and the IP-EOM Hamiltonian for the singles and
             doubles.
         """
+        if not ints:
+            ints = self.hbar_matvec_ip_intermediates(eris=eris, amplitudes=amplitudes)
         func, kwargs = self._load_function(
             "hbar_matvec_ip",
             eris=eris,
             amplitudes=amplitudes,
             excitations=excitations,
+            ints=ints,
         )
         res: Namespace[SpinArrayType] = func(**kwargs)
         res = util.Namespace(**{key.rstrip("new"): val for key, val in res.items()})
         return res
 
+    def hbar_matvec_ea_intermediates(
+        self,
+        eris: Optional[ERIsInputType] = None,
+        amplitudes: Optional[Namespace[SpinArrayType]] = None,
+    ) -> Namespace[NDArray[float]]:
+        """Compute intermediates for the EA-EOM Hamiltonian.
+
+        Args:
+            eris: Electron repulsion integrals.
+            amplitudes: Cluster amplitudes.
+
+        Returns:
+            Intermediate products for the EA-EOM Hamiltonian.
+        """
+        func, kwargs = self._load_function(
+            "hbar_matvec_ea_intermediates",
+            eris=eris,
+            amplitudes=amplitudes,
+        )
+        res: Namespace[NDArray[float]] = util.Namespace(**func(**kwargs))
+        return res
+
     def hbar_matvec_ea(
         self,
         excitations: Namespace[SpinArrayType],
         eris: Optional[ERIsInputType] = None,
         amplitudes: Optional[Namespace[SpinArrayType]] = None,
+        ints: Optional[Namespace[NDArray[float]]] = None,
     ) -> Namespace[SpinArrayType]:
         """Compute the product between a state vector and the EA-EOM Hamiltonian.
 
@@ -821,26 +871,53 @@ def hbar_matvec_ea(
             excitations: State vector as a set of excitation amplitudes.
             eris: Electron repulsion integrals.
             amplitudes: Cluster amplitudes.
+            ints: Intermediate products.
 
         Returns:
             Products between the state vectors and the EA-EOM Hamiltonian for the singles and
             doubles.
         """
+        if not ints:
+            ints = self.hbar_matvec_ea_intermediates(eris=eris, amplitudes=amplitudes)
         func, kwargs = self._load_function(
             "hbar_matvec_ea",
             eris=eris,
             amplitudes=amplitudes,
             excitations=excitations,
+            ints=ints,
         )
         res: Namespace[SpinArrayType] = func(**kwargs)
         res = util.Namespace(**{key.rstrip("new"): val for key, val in res.items()})
         return res
 
+    def hbar_matvec_ee_intermediates(
+        self,
+        eris: Optional[ERIsInputType] = None,
+        amplitudes: Optional[Namespace[SpinArrayType]] = None,
+    ) -> Namespace[NDArray[float]]:
+        """Compute intermediates for the EE-EOM Hamiltonian.
+
+        Args:
+            eris: Electron repulsion integrals.
+            amplitudes: Cluster amplitudes.
+
+        Returns:
+            Intermediate products for the EE-EOM Hamiltonian.
+        """
+        func, kwargs = self._load_function(
+            "hbar_matvec_ee_intermediates",
+            eris=eris,
+            amplitudes=amplitudes,
+        )
+        res: Namespace[NDArray[float]] = util.Namespace(**func(**kwargs))
+        return res
+
     def hbar_matvec_ee(
         self,
         excitations: Namespace[SpinArrayType],
         eris: Optional[ERIsInputType] = None,
         amplitudes: Optional[Namespace[SpinArrayType]] = None,
+        ints: Optional[Namespace[NDArray[float]]] = None,
     ) -> Namespace[SpinArrayType]:
         """Compute the product between a state vector and the EE-EOM Hamiltonian.
 
@@ -848,16 +925,20 @@ def hbar_matvec_ee(
             excitations: State vector as a set of excitation amplitudes.
             eris: Electron repulsion integrals.
             amplitudes: Cluster amplitudes.
+            ints: Intermediate products.
 
         Returns:
             Products between the state vectors and the EE-EOM Hamiltonian for the singles and
             doubles.
         """
+        if not ints:
+            ints = self.hbar_matvec_ee_intermediates(eris=eris, amplitudes=amplitudes)
         func, kwargs = self._load_function(
             "hbar_matvec_ee",
             eris=eris,
             amplitudes=amplitudes,
             excitations=excitations,
+            ints=ints,
         )
         res: Namespace[SpinArrayType] = func(**kwargs)
         res = util.Namespace(**{key.rstrip("new"): val for key, val in res.items()})
diff --git a/ebcc/eom/base.py b/ebcc/eom/base.py
index fe19ce04..46d55ebd 100644
--- a/ebcc/eom/base.py
+++ b/ebcc/eom/base.py
@@ -141,19 +141,37 @@ def vector_to_amplitudes(self, vector: NDArray[float]) -> Namespace[SpinArrayTyp
 
     @abstractmethod
     def matvec(
-        self, vector: NDArray[float], eris: Optional[ERIsInputType] = None
+        self,
+        vector: NDArray[float],
+        eris: Optional[ERIsInputType] = None,
+        ints: Optional[NDArray[float]] = None,
     ) -> NDArray[float]:
         """Apply the Hamiltonian to a vector.
 
         Args:
             vector: State vector to apply the Hamiltonian to.
             eris: Electronic repulsion integrals.
+            ints: Intermediate products.
 
         Returns:
             Resulting vector.
         """
         pass
 
+    @abstractmethod
+    def matvec_intermediates(
+        self, eris: Optional[ERIsInputType] = None
+    ) -> Namespace[NDArray[float]]:
+        """Get the intermediates for application of the Hamiltonian to a vector.
+
+        Args:
+            eris: Electronic repulsion integrals.
+
+        Returns:
+            Intermediate products.
+        """
+        pass
+
     @abstractmethod
     def diag(self, eris: Optional[ERIsInputType] = None) -> NDArray[float]:
         """Get the diagonal of the Hamiltonian.
@@ -263,7 +281,8 @@ def davidson(
         )
 
         # Get the matrix-vector products and the diagonal:
-        matvecs = lambda vs: [self.matvec(v, eris=eris) for v in vs]
+        ints = self.matvec_intermediates(eris=eris)
+        matvecs = lambda vs: [self.matvec(v, eris=eris, ints=ints) for v in vs]
         diag = self.diag(eris=eris)
 
         # Get the guesses:
@@ -367,21 +386,38 @@ def vector_to_amplitudes(self, vector: NDArray[float]) -> Namespace[SpinArrayTyp
         return self.ebcc.vector_to_excitations_ip(vector)
 
     def matvec(
-        self, vector: NDArray[float], eris: Optional[ERIsInputType] = None
+        self,
+        vector: NDArray[float],
+        eris: Optional[ERIsInputType] = None,
+        ints: Optional[NDArray[float]] = None,
     ) -> NDArray[float]:
         """Apply the Hamiltonian to a vector.
 
         Args:
             vector: State vector to apply the Hamiltonian to.
             eris: Electronic repulsion integrals.
+            ints: Intermediate products.
 
         Returns:
             Resulting vector.
         """
         amplitudes = self.vector_to_amplitudes(vector)
-        result = self.ebcc.hbar_matvec_ip(amplitudes, eris=eris)
+        result = self.ebcc.hbar_matvec_ip(amplitudes, eris=eris, ints=ints)
         return self.amplitudes_to_vector(result)
 
+    def matvec_intermediates(
+        self, eris: Optional[ERIsInputType] = None
+    ) -> Namespace[NDArray[float]]:
+        """Get the intermediates for application of the Hamiltonian to a vector.
+
+        Args:
+            eris: Electronic repulsion integrals.
+
+        Returns:
+            Intermediate products.
+        """
+        return self.ebcc.hbar_matvec_ip_intermediates(eris=eris)
+
 
 class BaseEA_EOM(BaseEOM):
     """Base class for electron-affinity EOM-CC."""
@@ -414,21 +450,38 @@ def vector_to_amplitudes(self, vector: NDArray[float]) -> Namespace[SpinArrayTyp
         return self.ebcc.vector_to_excitations_ea(vector)
 
     def matvec(
-        self, vector: NDArray[float], eris: Optional[ERIsInputType] = None
+        self,
+        vector: NDArray[float],
+        eris: Optional[ERIsInputType] = None,
+        ints: Optional[NDArray[float]] = None,
     ) -> NDArray[float]:
         """Apply the Hamiltonian to a vector.
 
         Args:
             vector: State vector to apply the Hamiltonian to.
             eris: Electronic repulsion integrals.
+            ints: Intermediate products.
 
         Returns:
             Resulting vector.
         """
         amplitudes = self.vector_to_amplitudes(vector)
-        result = self.ebcc.hbar_matvec_ea(amplitudes, eris=eris)
+        result = self.ebcc.hbar_matvec_ea(amplitudes, eris=eris, ints=ints)
         return self.amplitudes_to_vector(result)
 
+    def matvec_intermediates(
+        self, eris: Optional[ERIsInputType] = None
+    ) -> Namespace[NDArray[float]]:
+        """Get the intermediates for application of the Hamiltonian to a vector.
+
+        Args:
+            eris: Electronic repulsion integrals.
+
+        Returns:
+            Intermediate products.
+        """
+        return self.ebcc.hbar_matvec_ea_intermediates(eris=eris)
+
 
 class BaseEE_EOM(BaseEOM):
     """Base class for electron-electron EOM-CC."""
@@ -461,17 +514,34 @@ def vector_to_amplitudes(self, vector: NDArray[float]) -> Namespace[SpinArrayTyp
         return self.ebcc.vector_to_excitations_ee(vector)
 
     def matvec(
-        self, vector: NDArray[float], eris: Optional[ERIsInputType] = None
+        self,
+        vector: NDArray[float],
+        eris: Optional[ERIsInputType] = None,
+        ints: Optional[NDArray[float]] = None,
     ) -> NDArray[float]:
         """Apply the Hamiltonian to a vector.
 
         Args:
             vector: State vector to apply the Hamiltonian to.
             eris: Electronic repulsion integrals.
+            ints: Intermediate products.
 
         Returns:
             Resulting vector.
         """
         amplitudes = self.vector_to_amplitudes(vector)
-        result = self.ebcc.hbar_matvec_ee(amplitudes, eris=eris)
+        result = self.ebcc.hbar_matvec_ee(amplitudes, eris=eris, ints=ints)
         return self.amplitudes_to_vector(result)
+
+    def matvec_intermediates(
+        self, eris: Optional[ERIsInputType] = None
+    ) -> Namespace[NDArray[float]]:
+        """Get the intermediates for application of the Hamiltonian to a vector.
+
+        Args:
+            eris: Electronic repulsion integrals.
+
+        Returns:
+            Intermediate products.
+        """
+        return self.ebcc.hbar_matvec_ee_intermediates(eris=eris)
diff --git a/pyproject.toml b/pyproject.toml
index 134e5fc8..669f1837 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "ebcc"
-version = "1.4.5"
+version = "1.5.0"
 description = "Coupled cluster calculations on electron-boson systems"
 keywords = [
     "quantum", "chemistry",