diff --git a/.gitignore b/.gitignore
index 3557f10..2462954 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,6 @@
 __pycache__/
+.pytest_cache/
 .DS_STORE
 .vscode
-*.egg-info
\ No newline at end of file
+*.egg-info
+build/
\ No newline at end of file
diff --git a/adcgen/__init__.py b/adcgen/__init__.py
index 3343cae..f2aad38 100644
--- a/adcgen/__init__.py
+++ b/adcgen/__init__.py
@@ -20,6 +20,7 @@
 from .sympy_objects import (AntiSymmetricTensor, SymmetricTensor, Amplitude,
                             NonSymmetricTensor, KroneckerDelta, SymbolicTensor)
 from .tensor_names import tensor_names
+from .resolution_of_identity import apply_resolution_of_identity
 from . import sort_expr as sort
 
 
@@ -35,6 +36,7 @@
            "Intermediates", "reduce_expr", "factor_intermediates",
            "sort",
            "transform_to_spatial_orbitals",
+           "apply_resolution_of_identity",
            "apply_cvs_approximation",
            "generate_code", "optimize_contractions",
            "unoptimized_contraction", "Contraction",
diff --git a/adcgen/expression/expr_container.py b/adcgen/expression/expr_container.py
index 7a6e216..e28f10e 100644
--- a/adcgen/expression/expr_container.py
+++ b/adcgen/expression/expr_container.py
@@ -271,6 +271,29 @@ def rename_tensor(self, current: str, new: str) -> "ExprContainer":
         self._inner = renamed
         return self
 
+    def expand_coulomb_ri(self, factorisation: str = 'sym'
+                          ) -> 'ExprContainer':
+        """
+        Expands the Coulomb operators (pq | rs) into RI format
+
+        Parameters
+        ----------
+        factorisation : str, optional
+            The type of factorisation ('sym' or 'asym'), by default 'sym'
+
+        Returns
+        -------
+        ExprContainer
+            The factorised expression.
+        """
+        res = S.Zero
+        for term in self.terms:
+            res += term.expand_coulomb_ri(factorisation=factorisation,
+                                          wrap_result=False)
+        assert isinstance(res, Expr)
+        self._inner = res
+        return self
+
     def expand_antisym_eri(self) -> "ExprContainer":
         """
         Expands the antisymmetric ERI using chemists notation
diff --git a/adcgen/expression/object_container.py b/adcgen/expression/object_container.py
index 24f0f57..cb39958 100644
--- a/adcgen/expression/object_container.py
+++ b/adcgen/expression/object_container.py
@@ -445,6 +445,11 @@ def allowed_spin_blocks(self) -> tuple[str, ...] | None:
                 ]))
             elif name == tensor_names.coulomb:  # ERI in chemist notation
                 return ("aaaa", "aabb", "bbaa", "bbbb")
+            elif name in (tensor_names.ri_sym, tensor_names.ri_asym_eri,
+                          tensor_names.ri_asym_factor):
+                return ("aaa", "abb")
+            elif name == tensor_names.fock:
+                return ("aa", "bb")
         elif isinstance(obj, KroneckerDelta):  # delta
             # spins have to be equal
             return ("aa", "bb")
@@ -499,6 +504,8 @@ def format_indices(indices) -> str:
                 ),
                 # coulomb integral chemist notation
                 tensor_names.coulomb: lambda up, lo: f"({up}\\vert {lo})",
+                # 2e3c integral in asymmetric RI
+                tensor_names.ri_asym_eri: lambda up, lo: f"({up}\\vert {lo})",
                 # orbital energy
                 tensor_names.orb_energy: lambda _, lo: f"\\varepsilon_{{{lo}}}"
             }
@@ -724,6 +731,84 @@ def rename_tensor(self, current: str, new: str,
             obj = ExprContainer(obj, **self.assumptions)
         return obj
 
+    def expand_coulomb_ri(self, factorisation: str = 'sym',
+                          wrap_result: bool = True) -> "ExprContainer | Expr":
+        """
+        Expands the Coulomb operators (pq | rs) into RI format
+
+        Parameters
+        ----------
+        factorisation : str, optional
+            The type of factorisation ('sym' or 'asym'), by default 'sym'
+        wrap_result : bool, optional
+            Whether to wrap the result in an ExprContainer, by default True
+
+        Returns
+        -------
+        ExprContainer | Expr
+            The factorised expression.
+        """
+        from .expr_container import ExprContainer
+
+        if factorisation not in ('sym', 'asym'):
+            raise NotImplementedError("Only symmetric ('sym') and asymmetric "
+                                      "('asym') factorisation of the Coulomb "
+                                      "integral is implemented")
+
+        res = self.inner
+        base, exponent = self.base_and_exponent
+        if isinstance(base, SymmetricTensor) and \
+                base.name == tensor_names.coulomb:
+            if base.bra_ket_sym != 1:
+                raise NotImplementedError("Can only apply RI approximation to "
+                                          "coulomb integrals with "
+                                          "bra-ket symmetry.")
+            # we dont expand negative exponents, because the result
+            # (ab)^-n will evaluate to a^-n b^-n, which is
+            # only correct if the product ab has no contracted
+            # indices
+            if not exponent.is_Integer or exponent < S.Zero:
+                raise NotImplementedError("Can only apply RI approximation to "
+                                          "coulomb integrals "
+                                          "with positive integer exponents. "
+                                          f"{self} has an invalid exponent.")
+            # setup the assumptions for the aux index:
+            # assign alpha spin if represented in spatial orbitals
+            idx = self.idx
+            has_spin = bool(idx[0].spin)
+            if any(bool(s) != has_spin for s in idx):
+                raise NotImplementedError(f"The coulomb integral {self} has "
+                                          "to be represented either in spatial"
+                                          " or spin orbitals. A mixture is not"
+                                          " valid.")
+            assumptions = {"aux": True}
+            if has_spin:
+                assumptions["alpha"] = True
+            # actually do the expansion
+            p, q, r, s = idx
+            res = S.One
+            for _ in range(int(exponent)):  # exponent has to be positive
+                aux_idx = Index("P", **assumptions)
+                if factorisation == "sym":
+                    res *= SymmetricTensor(
+                        tensor_names.ri_sym, (aux_idx,), (p, q), 0
+                    )
+                    res *= SymmetricTensor(
+                        tensor_names.ri_sym, (aux_idx,), (r, s), 0
+                    )
+                else:
+                    assert factorisation == "asym"
+                    res *= SymmetricTensor(
+                        tensor_names.ri_asym_factor, (aux_idx,), (p, q), 0
+                    )
+                    res *= SymmetricTensor(
+                        tensor_names.ri_asym_eri, (aux_idx,), (r, s), 0
+                    )
+        if wrap_result:
+            kwargs = self.assumptions
+            res = ExprContainer(res, **kwargs)
+        return res
+
     def expand_antisym_eri(self, wrap_result: bool = True
                            ) -> "ExprContainer | Expr":
         """
@@ -782,6 +867,12 @@ def expand_intermediates(self, target: Sequence[Index],
             False: The intermediate is only expanded once, e.g., n'th
               order MP t-amplitudes are expressed by means of (n-1)'th order
               MP t-amplitudes and ERI.
+        braket_sym_tensors: Sequence[str], optional
+            Add bra-ket-symmetry to the given tensors of the expanded
+            expression (after expansion of the intermediates).
+        braket_antisym_tensors: Sequence[str], optional
+            Add bra-ket-antisymmetry to the given tensors of the expanded
+            expression (after expansion of the intermediates).
         """
         from ..intermediates import Intermediates
         from .expr_container import ExprContainer
@@ -800,18 +891,31 @@ def expand_intermediates(self, target: Sequence[Index],
         itmd = Intermediates().available.get(longname, None)
         expanded = self.inner
         if itmd is not None:
+            # for negative exponents we would have to ensure that the
+            # intermediate is a "long" intermediate that consists of
+            # multiple terms. Or if it consists of a single term
+            # that it does not have any contracted indices
+            # However, this can only be checked by calling ".expand()"
+            # on the contributions in the for loop below, which seems bad.
+            # A short intermediates will be expanded as
+            # X^-2 = (ab * cd)^-1 -> a^-1 b^-1 c^-1 d^-1
+            # where the last step is not correct if the intermediate
+            # has contracted indices.
+            exponent = self.exponent
+            if not exponent.is_Integer or exponent < S.Zero:
+                raise NotImplementedError(
+                    "Can only expand intermediates with positive "
+                    f"integer exponents. {self} has an invalid exponent."
+                )
             # Use a for loop to obtain different contracted itmd indices
             # for each x in: x * x * ...
             expanded = S.One
-            exponent = self.exponent
             assert exponent.is_Integer
             for _ in range(abs(int(exponent))):
                 expanded *= itmd.expand_itmd(
                     indices=self.idx, wrap_result=False,
                     fully_expand=fully_expand
                 )
-            if exponent < S.Zero:
-                expanded = Pow(expanded, -1)
         # apply assumptions to the expanded object
         if braket_sym_tensors or braket_antisym_tensors:
             expanded = ExprContainer(expanded).add_bra_ket_sym(
diff --git a/adcgen/expression/polynom_container.py b/adcgen/expression/polynom_container.py
index 57ab724..b841241 100644
--- a/adcgen/expression/polynom_container.py
+++ b/adcgen/expression/polynom_container.py
@@ -2,7 +2,7 @@
 from functools import cached_property
 from typing import Any, TYPE_CHECKING
 
-from sympy import Add, Expr, Pow, S
+from sympy import Add, Expr, Mul, Pow, S
 
 from ..indices import Index, sort_idx_canonical
 from .container import Container
@@ -188,6 +188,46 @@ def rename_tensor(self, current: str, new: str,
             renamed = ExprContainer(inner=renamed, **self.assumptions)
         return renamed
 
+    def expand_coulomb_ri(self, factorisation: str = 'sym',
+                          wrap_result: bool = True) -> "Expr | ExprContainer":
+        """
+        Expands the Coulomb operators (pq | rs) into RI format
+
+        Parameters
+        ----------
+        factorisation : str, optional
+            The type of factorisation ('sym' or 'asym'), by default 'sym'
+        wrap_result : bool, optional
+            Whether to wrap the result in an ExprContainer, by default True
+        """
+        from .expr_container import ExprContainer
+
+        exponent = self.exponent
+        if not exponent.is_Integer:
+            raise ValueError("Can only apply RI approximation to Polynomials "
+                             "with integer exponents. "
+                             f"{self} has an invalid exponent.")
+        # use a for loop so the contracted aux indices for each x in
+        # x * x * ... = x^n are different.
+        expanded = S.One
+        for _ in range(abs(int(exponent))):
+            contrib = S.Zero
+            for term in self.terms:
+                contrib += term.expand_coulomb_ri(
+                    factorisation=factorisation, wrap_result=False
+                )
+            assert isinstance(contrib, Expr)
+            if exponent < S.Zero:
+                # a mul object would be simplified as
+                # (ab)^-1 -> a^-1 b^-1
+                # which is only correct if a and b have no contracted indices.
+                assert not isinstance(contrib, Mul)
+                contrib = Pow(contrib, -1)
+            expanded *= contrib
+        if wrap_result:
+            expanded = ExprContainer(inner=expanded, **self.assumptions)
+        return expanded
+
     def expand_antisym_eri(self, wrap_result: bool = True):
         """
         Expands the antisymmetric ERI using chemists notation
@@ -217,16 +257,30 @@ def expand_intermediates(self, target: Sequence[Index],
         """Expands all known intermediates in the polynom."""
         from .expr_container import ExprContainer
 
-        expanded = S.Zero
-        for term in self.terms:
-            expanded += term.expand_intermediates(
-                target, wrap_result=False, fully_expand=fully_expand,
-                braket_sym_tensors=braket_sym_tensors,
-                braket_antisym_tensors=braket_antisym_tensors
-            )
-        assert isinstance(expanded, Expr)
-        expanded = Pow(expanded, self.exponent)
-
+        exponent = self.exponent
+        if not exponent.is_Integer:
+            raise NotImplementedError("Can only expand intermediates for "
+                                      "polynoms with integer exponents."
+                                      f"{self} has an invalid exponent.")
+        # use a for loop so the contracted itmd indices for each x in
+        # x * x * ... = x^n are different.
+        expanded = S.One
+        for _ in range(abs(int(exponent))):
+            contrib = S.Zero
+            for term in self.terms:
+                contrib += term.expand_intermediates(
+                    target, wrap_result=False, fully_expand=fully_expand,
+                    braket_sym_tensors=braket_sym_tensors,
+                    braket_antisym_tensors=braket_antisym_tensors
+                )
+            assert isinstance(contrib, Expr)
+            if exponent < S.Zero:
+                # a mul object would be simplified as
+                # (ab)^-1 -> a^-1 b^-1
+                # which is only correct if a and b have no contracted indices.
+                assert not isinstance(contrib, Mul)
+                contrib = Pow(contrib, -1)
+            expanded *= contrib
         if wrap_result:
             assumptions = self.assumptions
             assumptions["target_idx"] = target
diff --git a/adcgen/expression/term_container.py b/adcgen/expression/term_container.py
index c23e060..348bf31 100644
--- a/adcgen/expression/term_container.py
+++ b/adcgen/expression/term_container.py
@@ -583,6 +583,35 @@ def rename_tensor(self, current: str, new: str, wrap_result: bool = True
             renamed = ExprContainer(renamed, **self.assumptions)
         return renamed
 
+    def expand_coulomb_ri(self, factorisation: str = 'sym',
+                          wrap_result: bool = True) -> "Expr | ExprContainer":
+        """
+        Expands the Coulomb operators (pq | rs) into RI format
+
+        Parameters
+        ----------
+        factorisation : str, optional
+            The type of factorisation ('sym' or 'asym'), by default 'sym'
+        wrap_result : bool, optional
+            Whether to wrap the result in an ExprContainer, by default True
+
+        Returns
+        -------
+        ExprContainer | Expr
+            The factorised expression.
+        """
+        from .expr_container import ExprContainer
+
+        factorised = S.One
+        for obj in self.objects:
+            factorised *= obj.expand_coulomb_ri(factorisation=factorisation,
+                                                wrap_result=False)
+
+        if wrap_result:
+            assumptions = self.assumptions
+            factorised = ExprContainer(factorised, **assumptions)
+        return factorised
+
     def expand_antisym_eri(self, wrap_result: bool = True
                            ) -> "ExprContainer | Expr":
         """
diff --git a/adcgen/func.py b/adcgen/func.py
index ae8e5ad..ec5319b 100644
--- a/adcgen/func.py
+++ b/adcgen/func.py
@@ -159,7 +159,10 @@ def import_tensor(tensor: str) -> Expr:
             # ADC-Amplitude or t-amplitudes
             if is_adc_amplitude(name) or is_t_amplitude(name):
                 base: Expr = Amplitude(name, upper, lower)
-            elif name == tensor_names.coulomb:  # eri in chemist notation
+            elif name in (tensor_names.coulomb, tensor_names.ri_sym,
+                          tensor_names.ri_asym_eri,
+                          tensor_names.ri_asym_factor):
+                # eri in chemist notation or RI tensor
                 base: Expr = SymmetricTensor(name, upper, lower)
             else:
                 base: Expr = AntiSymmetricTensor(name, upper, lower)
diff --git a/adcgen/generate_code/config.json b/adcgen/generate_code/config.json
index 676fbe8..0876473 100644
--- a/adcgen/generate_code/config.json
+++ b/adcgen/generate_code/config.json
@@ -2,6 +2,7 @@
   "sizes": {
     "core": 5,
     "occ": 20,
-    "virt": 200
+    "virt": 200,
+    "aux": 250
   }
 }
\ No newline at end of file
diff --git a/adcgen/generate_code/contraction.py b/adcgen/generate_code/contraction.py
index df9b957..4579518 100644
--- a/adcgen/generate_code/contraction.py
+++ b/adcgen/generate_code/contraction.py
@@ -23,6 +23,7 @@ class Sizes:
     occ: int = 0
     virt: int = 0
     general: int = 0
+    aux: int = 0
 
     @staticmethod
     def from_dict(input: dict[str, int]) -> "Sizes":
@@ -32,7 +33,7 @@ def from_dict(input: dict[str, int]) -> "Sizes":
         if not provided.
         """
         if "general" not in input:
-            input["general"] = sum(input.values())
+            input["general"] = sum(v for k, v in input.items() if k != "aux")
         return Sizes(**input)
 
     @staticmethod
@@ -232,6 +233,7 @@ class ScalingComponent:
     virt: int
     occ: int
     core: int
+    aux: int
 
     def evaluate_costs(self, sizes: Sizes) -> int:
         """
diff --git a/adcgen/indices.py b/adcgen/indices.py
index a281292..ab61638 100644
--- a/adcgen/indices.py
+++ b/adcgen/indices.py
@@ -45,6 +45,8 @@ def space(self) -> str:
             return "virt"
         elif self.assumptions0.get("core"):
             return "core"
+        elif self.assumptions0.get("aux"):
+            return "aux"
         else:
             return "general"
 
@@ -83,7 +85,7 @@ class Indices(metaclass=Singleton):
     # the valid spaces with their corresponding associated index names
     base = {
         "occ": "ijklmno", "virt": "abcdefgh", "general": "pqrstuvw",
-        "core": "IJKLMNO"
+        "core": "IJKLMNO", "aux": "PQRST"
     }
     # the valid spins
     spins = ("", "a", "b")
@@ -244,6 +246,8 @@ def _new_symbol(self, name: str, space: str, spin: str) -> Index:
             assumptions["above_fermi"] = True
         elif space == "core":
             assumptions["core"] = True
+        elif space == "aux":
+            assumptions["aux"] = True
         elif space != "general":
             raise ValueError(f"Invalid space {space}")
         if spin:
@@ -270,7 +274,7 @@ def sort_idx_canonical(idx: Index | Any):
         # - also add the hash here for wicks, where multiple i are around
         # - we have to map the spaces onto numbers, since in adcman and adcc
         # the ordering o < c < v is used for the definition of canonical blocks
-        space_keys = {"g": 0, "o": 1, "c": 2, "v": 3}
+        space_keys = {"g": 0, "o": 1, "c": 2, "v": 3, "a": 4}
         return (space_keys[idx.space[0]],
                 idx.spin,
                 int(idx.name[1:]) if idx.name[1:] else 0,
diff --git a/adcgen/resolution_of_identity.py b/adcgen/resolution_of_identity.py
new file mode 100644
index 0000000..d38bb25
--- /dev/null
+++ b/adcgen/resolution_of_identity.py
@@ -0,0 +1,75 @@
+from .expression import ExprContainer
+from .tensor_names import tensor_names
+from .misc import Inputerror
+from sympy import Symbol
+
+
+def apply_resolution_of_identity(expr: ExprContainer,
+                                 factorisation: str = 'sym',
+                                 resolve_indices: bool = True
+                                 ) -> ExprContainer:
+    """
+    Applies the Resolution of Identity approximation (RI, sometimes also
+    called density fitting, DF) to an expression. This implies that every
+    Coulomb operator is replaced by its factorised form. Two types of
+    factorisation are supported: symmetric and asymmetric.
+    In the symmetric decomposition, a Coulomb operator is approximated as:
+
+    (pq | rs) ~ B^P_{pq} B^P_{rs}
+    B^P_{pq} = (P | Q)^{-1/2} (Q | pq)
+
+    This decomposition is the default. In the asymmetric factorisation, the
+    same Coulomb operator is approximated as:
+
+    (pq | rs) ~ C^P_{pq} (P | rs)
+    C^P_{pq} = (P | Q)^{-1} (Q | pq)
+
+    Note that the RI approximation is only meaningful on Coulomb operator.
+    Therefore, this routine will crash and exit if the given expression has
+    not been expanded before. All RI indices receive an alpha spin
+    by default if the expression has been spin-integrated and no spin
+    otherwise.
+
+    Parameters
+    ----------
+    expr : ExprContainer
+        The expression to be factorised into RI format.
+    factorisation : str, optional
+        Which type of factorisation to use.
+        If 'sym', the symmetric factorisation variant is employed.
+        If 'asym', the asymmetric factorisation variant is employed
+        instead, by default 'sym'
+    resolve_indices: bool, optional
+        Whether the indices should be resolved to unique indices
+        after applying the RI approximation. This is equivalent
+        to calling 'substitute_contracted()' afterwards.
+
+    Returns
+    -------
+    ExprContainer
+        The factorised expression
+
+    Raises
+    ------
+    Inputerror
+        If a factorisation other than 'sym' or 'asym' is provided
+    Inputerror
+        If the expression still contains antisymmetric ERIs.
+    """
+
+    # Check if a valid factorisation is given
+    if factorisation not in ('sym', 'asym'):
+        raise Inputerror('Only symmetric (sym) and asymmetric (asym) '
+                         'factorisation modes are supported. '
+                         f'Received: {factorisation}')
+
+    # Check whether the expression contains antisymmetric ERIs
+    if Symbol(tensor_names.eri) in expr.inner.atoms(Symbol):
+        raise Inputerror('Resolution of Identity requires that the ERIs'
+                         ' be expanded first.')
+
+    ri_expr = expr.expand_coulomb_ri(factorisation=factorisation)
+    if resolve_indices:
+        ri_expr.substitute_contracted()
+
+    return ri_expr
diff --git a/adcgen/simplify.py b/adcgen/simplify.py
index 3ac02e7..9196858 100644
--- a/adcgen/simplify.py
+++ b/adcgen/simplify.py
@@ -60,7 +60,7 @@ def check_term(term: TermContainer) -> bool:
         # True if all requested tensors are in the term
         if strict == 'low':
             return all(t in available for t in set(t_strings))
-        # True if all requested Tensors occure the correct amount of times
+        # True if all requested Tensors occur the correct amount of times
         elif strict == 'medium':
             available = Counter(available)
             desired = Counter(t_strings)
diff --git a/adcgen/spatial_orbitals.py b/adcgen/spatial_orbitals.py
index 732df87..62881e3 100644
--- a/adcgen/spatial_orbitals.py
+++ b/adcgen/spatial_orbitals.py
@@ -179,11 +179,19 @@ def integrate_spin(expr: ExprContainer, target_idx: str,
                 addition: list[str | None] = [
                     None for _ in range(len(term_indices))
                 ]
+                # We keep an explicit note whether to keep the block
+                # This is required because of the internal loop structure
+                # of the checks
+                accept_spin_block: bool = True
                 for spin, idx in zip(block, indices):
                     if addition[idx] is not None and addition[idx] != spin:
-                        raise ValueError("Found invalid allowed spin block "
-                                         f"{block} for {obj}.")
+                        # This can occur if the same index appears in the same
+                        # object multiple times, e. g. < ij || ij >
+                        accept_spin_block = False
+                        break
                     addition[idx] = spin
+                if not accept_spin_block:
+                    continue
                 # check for contracdictions with the target_spin and skip the
                 # block if this is the case
                 if any(sp1 != sp2 for sp1, sp2 in
diff --git a/adcgen/sympy_objects.py b/adcgen/sympy_objects.py
index 6319218..9374612 100644
--- a/adcgen/sympy_objects.py
+++ b/adcgen/sympy_objects.py
@@ -320,9 +320,12 @@ def eval(cls, i: Expr, j: Expr) -> Expr | None:  # type: ignore[override]
             return S.Zero
 
         spi, spj = i.space[0], j.space[0]
-        valid_spaces = ["o", "v", "g", "c"]
+        valid_spaces = ["o", "v", "g", "c", "a"]
         assert spi in valid_spaces and spj in valid_spaces
-        if spi != "g" and spj != "g" and spi != spj:  # delta_ov / delta_vo
+        # delta_ov / delta_vo / delta_{aux general} / delta_{general aux}
+        if (spi != "g" and spj != "g" and spi != spj) or \
+                (spi == "g" and spj == "a") or \
+                (spi == "a" and spj == "g"):
             return S.Zero
         spi, spj = i.spin, j.spin
         assert spi in ["", "a", "b"] and spj in ["", "a", "b"]
@@ -367,12 +370,14 @@ def preferred_and_killable(self) -> tuple[Index, Index] | None:
         space2, spin2 = j.space[0], j.spin
         # ensure we have no unexpected space and spin
         assert (
-            space1 in ["o", "v", "g", "c"] and space2 in ["o", "v", "g", "c"]
+            space1 in ["o", "v", "g", "c", "a"]
+            and space2 in ["o", "v", "g", "c", "a"]
         )
         assert spin1 in ["", "a", "b"] and spin2 in ["", "a", "b"]
 
         if spin1 == spin2:  # nn / aa / bb  -> equal spin information
-            # oo / vv / cc / gg / og / vg / cg
+            # oo / vv / cc / gg / og / vg / cg / aa
+            # RI indices will always end up here
             if space1 == space2 or space2 == "g":
                 return (i, j)
             else:  # go / gv / gc
diff --git a/adcgen/tensor_names.json b/adcgen/tensor_names.json
index 1f68b32..38dd928 100644
--- a/adcgen/tensor_names.json
+++ b/adcgen/tensor_names.json
@@ -8,5 +8,8 @@
   "left_adc_amplitude": "X",
   "right_adc_amplitude": "Y",
   "orb_energy": "e",
-  "sym_orb_denom": "D"
+  "sym_orb_denom": "D",
+  "ri_sym": "B",
+  "ri_asym_factor": "C",
+  "ri_asym_eri": "G"
 }
\ No newline at end of file
diff --git a/adcgen/tensor_names.py b/adcgen/tensor_names.py
index e795724..b04ed3e 100644
--- a/adcgen/tensor_names.py
+++ b/adcgen/tensor_names.py
@@ -29,6 +29,15 @@ class TensorNames(metaclass=Singleton):
     the attributes storing the names are given in brackets:
     - antisymmetric ERI in physicist notation (eri): V
     - Coulomb integrals in chemist notation (coulomb): v
+    - Symmetrically decomposed RI integrals (ri_sym): B
+      These are formally calculated by decomposing the symmetric
+      four center integrals: (pq | rs) = B^Q_pq B^Q_rs
+    - The "factor" for an asymmetric RI integral (ri_asym_factor): C
+      This tensor is the dimensionless part of the asymmetric
+      resolution of identity decomposition: C^P_{pq} = (pq | Q) (Q | P)^{-1}
+    - The pure 2e3c RI integral (ri_asym_eri): G
+      This tensor is combined with C^P_{pq} to reform the symmetric
+      four center integrals: C^P_{pq} G^P_{rs} = (pq | rs)
     - The fock matrix (fock): f
     - The arbitrary N-particle operator matrix (operator): d
     - Ground state amplitudes (gs_amplitude): t
@@ -47,6 +56,9 @@ class TensorNames(metaclass=Singleton):
     """
     eri: str = "V"
     coulomb: str = "v"
+    ri_sym: str = "B"
+    ri_asym_factor: str = "C"
+    ri_asym_eri: str = "G"
     fock: str = "f"
     operator: str = "d"
     gs_amplitude: str = "t"
diff --git a/examples/resolution_of_identity.py b/examples/resolution_of_identity.py
new file mode 100644
index 0000000..7460bd4
--- /dev/null
+++ b/examples/resolution_of_identity.py
@@ -0,0 +1,41 @@
+from adcgen import (
+    Operators, GroundState, ExprContainer, transform_to_spatial_orbitals,
+    apply_resolution_of_identity
+)
+
+# We first declare the Hamiltonian operator
+op = Operators()
+
+# We exemplify this at the MP2 and MP3 energies
+# For this, we first define the ground state
+gs = GroundState(op)
+
+# Next, we can calculate the MP2 and MP3 energy
+energy_mp2 = ExprContainer(gs.energy(2))
+energy_mp3 = ExprContainer(gs.energy(3))
+
+# RI is only valid for real orbitals, wherefore we have to make these
+# expressions real
+energy_mp2.make_real()
+energy_mp3.make_real()
+
+# These can now be spin-integrated
+energy_mp2 = transform_to_spatial_orbitals(energy_mp2, '', '',
+                                           restricted=False)
+energy_mp3 = transform_to_spatial_orbitals(energy_mp3, '', '',
+                                           restricted=False)
+
+# Lastly, we can apply the resolution of identity approximation
+# We will decompose the MP2 energy symmetrically:
+energy_mp2 = apply_resolution_of_identity(energy_mp2, factorisation='sym')
+# And the MP3 energy asymetrically:
+energy_mp3 = apply_resolution_of_identity(energy_mp3, factorisation='asym')
+
+# We can now print the result
+print("RI-MP2 Energy:\n")
+print(energy_mp2.to_latex_str(spin_as_overbar=True))
+print()
+
+print("RI-MP3 Energy:\n")
+print(energy_mp3.to_latex_str(spin_as_overbar=True))
+print()
diff --git a/tests/contraction_test.py b/tests/contraction_test.py
index f67f8b1..b48aaca 100644
--- a/tests/contraction_test.py
+++ b/tests/contraction_test.py
@@ -42,15 +42,15 @@ def test_scaling(self):
         target_indices = (i, j)
         contr = Contraction(indices, names, target_indices)
         scaling = contr.scaling
-        comp = ScalingComponent(3, 0, 0, 3, 0)
-        mem = ScalingComponent(2, 0, 0, 2, 0)
+        comp = ScalingComponent(3, 0, 0, 3, 0, 0)
+        mem = ScalingComponent(2, 0, 0, 2, 0, 0)
         assert scaling.computational == comp
         assert scaling.memory == mem
         assert scaling == Scaling(comp, mem)
 
     def test_sizes(self):
         # test the automatic evaluation of the size of the general space
-        sizes = {"occ": 1, "virt": 2, "core": 3}
+        sizes = {"occ": 1, "virt": 2, "core": 3, "aux": 0}
         res = Sizes.from_dict(sizes)
         sizes["general"] = 6
         assert sizes == asdict(res)
@@ -59,17 +59,17 @@ def test_sizes(self):
         assert sizes == asdict(res)
 
     def test_evalute_costs(self):
-        sizes = {"occ": 3, "virt": 5, "core": 2, "general": 7}
+        sizes = {"occ": 3, "virt": 5, "core": 2, "general": 7, "aux": 8}
         sizes = Sizes.from_dict(sizes)
-        comp = ScalingComponent(42, 1, 2, 3, 4)
-        mem = ScalingComponent(42, 4, 3, 2, 1)
+        comp = ScalingComponent(42, 1, 2, 3, 4, 5)
+        mem = ScalingComponent(42, 5, 4, 3, 2, 1)
         scaling = Scaling(comp, mem)
-        assert comp.evaluate_costs(sizes) == 75600
-        assert mem.evaluate_costs(sizes) == 5402250
-        assert scaling.evaluate_costs(sizes) == (75600, 5402250)
+        assert comp.evaluate_costs(sizes) == 2477260800
+        assert mem.evaluate_costs(sizes) == 9075780000
+        assert scaling.evaluate_costs(sizes) == (2477260800, 9075780000)
         # ensure that zero sized spaces are ignored
-        sizes = {"occ": 3, "virt": 5, "core": 0}  # general == 8
+        sizes = {"occ": 3, "virt": 5, "core": 0, "aux": 0}  # general == 8
         sizes = Sizes.from_dict(sizes)
         assert comp.evaluate_costs(sizes) == 5400
-        comp = ScalingComponent(42, 0, 1, 2, 3)
+        comp = ScalingComponent(42, 0, 1, 2, 3, 0)
         assert comp.evaluate_costs(sizes) == 45
diff --git a/tests/indices_test.py b/tests/indices_test.py
index 40827e0..dba0d84 100644
--- a/tests/indices_test.py
+++ b/tests/indices_test.py
@@ -8,10 +8,15 @@ def test_get_indices(self):
         assert idx.get_indices("ijk") == idx.get_indices("ijk")
         assert idx.get_indices("ijk", "aba") == idx.get_indices("ijk", "aba")
         assert idx.get_indices("a", "a") != idx.get_indices("a", "b")
+        assert idx.get_indices("PQ") == idx.get_indices("PQ")
         res = idx.get_indices("Ij", "ba")
         I, j = res[("core", "b")].pop(), res[("occ", "a")].pop()
         assert I.space == "core" and I.spin == "b"
         assert j.space == "occ" and j.spin == "a"
+        res = idx.get_indices("Pa")
+        P, a = res[("aux", "")].pop(), res[("virt", "")].pop()
+        assert P.space == "aux" and P.spin == ""
+        assert a.space == "virt" and a.spin == ""
 
     def test_get_generic_indices(self):
         # ensure that generic indices don't overlap
diff --git a/tests/reference_data/generate_data.py b/tests/reference_data/generate_data.py
index 7f5bdae..d89bdf5 100644
--- a/tests/reference_data/generate_data.py
+++ b/tests/reference_data/generate_data.py
@@ -10,6 +10,8 @@
 from adcgen.secular_matrix import SecularMatrix
 from adcgen.simplify import simplify, remove_tensor
 from adcgen.tensor_names import tensor_names
+from adcgen.resolution_of_identity import apply_resolution_of_identity
+from adcgen.spatial_orbitals import transform_to_spatial_orbitals
 from adcgen import sort_expr as sort
 
 import itertools
@@ -364,6 +366,48 @@ def gen_adc_properties_trans_moment(self):
                         dump["real_transition_dm"][block] = str(expr)
         write_json(results, outfile)
 
+    def gen_ri_gs_energy(self):
+        results: dict = {}
+        outfile = "ri_gs_energy.json"
+
+        variations = itertools.product(['mp', 're'], [0, 1, 2, 3], ['r', 'u'],
+                                       ['sym', 'asym'])
+
+        for variant, order, restriction, symmetry in variations:
+            if variant not in results:
+                results[variant] = {}
+            if order not in results[variant]:
+                results[variant][order] = {}
+            if restriction not in results[variant][order]:
+                results[variant][order][restriction] = {}
+            gs = self.gs[variant]
+            gs_energy = ExprContainer(gs.energy(order), real=True)
+            restricted = restriction == 'r'
+            gs_energy = transform_to_spatial_orbitals(gs_energy, '', '',
+                                                      restricted=restricted)
+            gs_energy = apply_resolution_of_identity(gs_energy, symmetry)
+            gs_energy.substitute_contracted()
+            results[variant][order][restriction][symmetry] = str(gs_energy)
+        write_json(results, outfile)
+
+    def gen_spatial_gs_energy(self):
+        outfile = "spatial_gs_energy.json"
+
+        results: dict = {}
+        for variant in ['mp', 're']:
+            results[variant] = {}
+            gs = self.gs[variant]
+            for order in [0, 1, 2, 3]:
+                results[variant][order] = {}
+                for restriction in ['r', 'u']:
+                    energy = ExprContainer(gs.energy(order), real=True)
+                    restr = restriction == 'r'
+                    energy = transform_to_spatial_orbitals(energy, '', '',
+                                                           restricted=restr)
+                    energy.substitute_contracted()
+                    results[variant][order][restriction] = str(energy)
+        write_json(results, outfile)
+
 
 def write_json(data, filename):
     json.dump(data, open(filename, 'w'), indent=2)
diff --git a/tests/reference_data/ri_gs_energy.json b/tests/reference_data/ri_gs_energy.json
new file mode 100644
index 0000000..bfa2d61
--- /dev/null
+++ b/tests/reference_data/ri_gs_energy.json
@@ -0,0 +1,86 @@
+{
+  "mp": {
+    "0": {
+      "r": {
+        "sym": "2 {f^{i_{\\alpha}}_{i_{\\alpha}}}",
+        "asym": "2 {f^{i_{\\alpha}}_{i_{\\alpha}}}"
+      },
+      "u": {
+        "sym": "{f^{i_{\\alpha}}_{i_{\\alpha}}} + {f^{i_{\\beta}}_{i_{\\beta}}}",
+        "asym": "{f^{i_{\\alpha}}_{i_{\\alpha}}} + {f^{i_{\\beta}}_{i_{\\beta}}}"
+      }
+    },
+    "1": {
+      "r": {
+        "sym": "- 2 {B^{P_{\\alpha}}_{i_{\\alpha}i_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}j_{\\alpha}}} + {B^{P_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}}^{2}",
+        "asym": "- 2 {C^{P_{\\alpha}}_{i_{\\alpha}i_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}j_{\\alpha}}} + {C^{P_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {G^{P_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}}"
+      },
+      "u": {
+        "sym": "- {B^{P_{\\alpha}}_{i_{\\alpha}i_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\beta}i_{\\beta}}} - \\frac{{B^{P_{\\alpha}}_{i_{\\alpha}i_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}j_{\\alpha}}}}{2} + \\frac{{B^{P_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}}^{2}}{2} - \\frac{{B^{P_{\\alpha}}_{i_{\\beta}i_{\\beta}}} {B^{P_{\\alpha}}_{j_{\\beta}j_{\\beta}}}}{2} + \\frac{{B^{P_{\\alpha}}_{i_{\\beta}j_{\\beta}}}^{2}}{2}",
+        "asym": "- {C^{P_{\\alpha}}_{i_{\\alpha}i_{\\alpha}}} {G^{P_{\\alpha}}_{i_{\\beta}i_{\\beta}}} - \\frac{{C^{P_{\\alpha}}_{i_{\\alpha}i_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}j_{\\alpha}}}}{2} + \\frac{{C^{P_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {G^{P_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}}}{2} - \\frac{{C^{P_{\\alpha}}_{i_{\\beta}i_{\\beta}}} {G^{P_{\\alpha}}_{j_{\\beta}j_{\\beta}}}}{2} + \\frac{{C^{P_{\\alpha}}_{i_{\\beta}j_{\\beta}}} {G^{P_{\\alpha}}_{i_{\\beta}j_{\\beta}}}}{2}"
+      }
+    },
+    "2": {
+      "r": {
+        "sym": "- \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{2} + \\frac{3 {t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{2}",
+        "asym": "- \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{2} + \\frac{3 {t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{2}"
+      },
+      "u": {
+        "sym": "- {t1^{a_{\\alpha}a_{\\beta}}_{i_{\\alpha}i_{\\beta}}} {B^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} - \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{4} + \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{4} - \\frac{{t1^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {B^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} {B^{P_{\\alpha}}_{j_{\\beta}b_{\\beta}}}}{4} + \\frac{{t1^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {B^{P_{\\alpha}}_{i_{\\beta}b_{\\beta}}} {B^{P_{\\alpha}}_{j_{\\beta}a_{\\beta}}}}{4}",
+        "asym": "- {t1^{a_{\\alpha}a_{\\beta}}_{i_{\\alpha}i_{\\beta}}} {C^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {G^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} - \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{4} + \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{4} - \\frac{{t1^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {C^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} {G^{P_{\\alpha}}_{j_{\\beta}b_{\\beta}}}}{4} + \\frac{{t1^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {C^{P_{\\alpha}}_{i_{\\beta}b_{\\beta}}} {G^{P_{\\alpha}}_{j_{\\beta}a_{\\beta}}}}{4}"
+      }
+    },
+    "3": {
+      "r": {
+        "sym": "- \\frac{3 {t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{2} + \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{2}",
+        "asym": "- \\frac{3 {t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{2} + \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{2}"
+      },
+      "u": {
+        "sym": "- {t2^{a_{\\alpha}a_{\\beta}}_{i_{\\alpha}i_{\\beta}}} {B^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} - \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{4} + \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{4} - \\frac{{t2^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {B^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} {B^{P_{\\alpha}}_{j_{\\beta}b_{\\beta}}}}{4} + \\frac{{t2^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {B^{P_{\\alpha}}_{i_{\\beta}b_{\\beta}}} {B^{P_{\\alpha}}_{j_{\\beta}a_{\\beta}}}}{4}",
+        "asym": "- {t2^{a_{\\alpha}a_{\\beta}}_{i_{\\alpha}i_{\\beta}}} {C^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {G^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} - \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{4} + \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{4} - \\frac{{t2^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {C^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} {G^{P_{\\alpha}}_{j_{\\beta}b_{\\beta}}}}{4} + \\frac{{t2^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {C^{P_{\\alpha}}_{i_{\\beta}b_{\\beta}}} {G^{P_{\\alpha}}_{j_{\\beta}a_{\\beta}}}}{4}"
+      }
+    }
+  },
+  "re": {
+    "0": {
+      "r": {
+        "sym": "2 {f^{i_{\\alpha}}_{i_{\\alpha}}} - 2 {B^{P_{\\alpha}}_{i_{\\alpha}i_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}j_{\\alpha}}} + {B^{P_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}}^{2}",
+        "asym": "2 {f^{i_{\\alpha}}_{i_{\\alpha}}} - 2 {C^{P_{\\alpha}}_{i_{\\alpha}i_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}j_{\\alpha}}} + {C^{P_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {G^{P_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}}"
+      },
+      "u": {
+        "sym": "{f^{i_{\\alpha}}_{i_{\\alpha}}} + {f^{i_{\\beta}}_{i_{\\beta}}} - {B^{P_{\\alpha}}_{i_{\\alpha}i_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\beta}i_{\\beta}}} - \\frac{{B^{P_{\\alpha}}_{i_{\\alpha}i_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}j_{\\alpha}}}}{2} + \\frac{{B^{P_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}}^{2}}{2} - \\frac{{B^{P_{\\alpha}}_{i_{\\beta}i_{\\beta}}} {B^{P_{\\alpha}}_{j_{\\beta}j_{\\beta}}}}{2} + \\frac{{B^{P_{\\alpha}}_{i_{\\beta}j_{\\beta}}}^{2}}{2}",
+        "asym": "{f^{i_{\\alpha}}_{i_{\\alpha}}} + {f^{i_{\\beta}}_{i_{\\beta}}} - {C^{P_{\\alpha}}_{i_{\\alpha}i_{\\alpha}}} {G^{P_{\\alpha}}_{i_{\\beta}i_{\\beta}}} - \\frac{{C^{P_{\\alpha}}_{i_{\\alpha}i_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}j_{\\alpha}}}}{2} + \\frac{{C^{P_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {G^{P_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}}}{2} - \\frac{{C^{P_{\\alpha}}_{i_{\\beta}i_{\\beta}}} {G^{P_{\\alpha}}_{j_{\\beta}j_{\\beta}}}}{2} + \\frac{{C^{P_{\\alpha}}_{i_{\\beta}j_{\\beta}}} {G^{P_{\\alpha}}_{i_{\\beta}j_{\\beta}}}}{2}"
+      }
+    },
+    "1": {
+      "r": {
+        "sym": "0",
+        "asym": "0"
+      },
+      "u": {
+        "sym": "0",
+        "asym": "0"
+      }
+    },
+    "2": {
+      "r": {
+        "sym": "- \\frac{3 {t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{2} + \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{2}",
+        "asym": "- \\frac{3 {t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{2} + \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{2}"
+      },
+      "u": {
+        "sym": "- {t1^{a_{\\alpha}a_{\\beta}}_{i_{\\alpha}i_{\\beta}}} {B^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} - \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{4} + \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{4} - \\frac{{t1^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {B^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} {B^{P_{\\alpha}}_{j_{\\beta}b_{\\beta}}}}{4} + \\frac{{t1^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {B^{P_{\\alpha}}_{i_{\\beta}b_{\\beta}}} {B^{P_{\\alpha}}_{j_{\\beta}a_{\\beta}}}}{4}",
+        "asym": "- {t1^{a_{\\alpha}a_{\\beta}}_{i_{\\alpha}i_{\\beta}}} {C^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {G^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} - \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{4} + \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{4} - \\frac{{t1^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {C^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} {G^{P_{\\alpha}}_{j_{\\beta}b_{\\beta}}}}{4} + \\frac{{t1^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {C^{P_{\\alpha}}_{i_{\\beta}b_{\\beta}}} {G^{P_{\\alpha}}_{j_{\\beta}a_{\\beta}}}}{4}"
+      }
+    },
+    "3": {
+      "r": {
+        "sym": "2 {t2^{a_{\\alpha}}_{i_{\\alpha}}} {f^{i_{\\alpha}}_{a_{\\alpha}}} - \\frac{3 {t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{2} + \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{2}",
+        "asym": "2 {t2^{a_{\\alpha}}_{i_{\\alpha}}} {f^{i_{\\alpha}}_{a_{\\alpha}}} - \\frac{3 {t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{2} + \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{2}"
+      },
+      "u": {
+        "sym": "{t2^{a_{\\alpha}}_{i_{\\alpha}}} {f^{i_{\\alpha}}_{a_{\\alpha}}} + {t2^{a_{\\beta}}_{i_{\\beta}}} {f^{i_{\\beta}}_{a_{\\beta}}} - {t2^{a_{\\alpha}a_{\\beta}}_{i_{\\alpha}i_{\\beta}}} {B^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} - \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{4} + \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {B^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {B^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{4} - \\frac{{t2^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {B^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} {B^{P_{\\alpha}}_{j_{\\beta}b_{\\beta}}}}{4} + \\frac{{t2^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {B^{P_{\\alpha}}_{i_{\\beta}b_{\\beta}}} {B^{P_{\\alpha}}_{j_{\\beta}a_{\\beta}}}}{4}",
+        "asym": "{t2^{a_{\\alpha}}_{i_{\\alpha}}} {f^{i_{\\alpha}}_{a_{\\alpha}}} + {t2^{a_{\\beta}}_{i_{\\beta}}} {f^{i_{\\beta}}_{a_{\\beta}}} - {t2^{a_{\\alpha}a_{\\beta}}_{i_{\\alpha}i_{\\beta}}} {C^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {G^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} - \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}a_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{4} + \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {C^{P_{\\alpha}}_{i_{\\alpha}b_{\\alpha}}} {G^{P_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{4} - \\frac{{t2^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {C^{P_{\\alpha}}_{i_{\\beta}a_{\\beta}}} {G^{P_{\\alpha}}_{j_{\\beta}b_{\\beta}}}}{4} + \\frac{{t2^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {C^{P_{\\alpha}}_{i_{\\beta}b_{\\beta}}} {G^{P_{\\alpha}}_{j_{\\beta}a_{\\beta}}}}{4}"
+      }
+    }
+  }
+}
\ No newline at end of file
diff --git a/tests/reference_data/spatial_gs_energy.json b/tests/reference_data/spatial_gs_energy.json
new file mode 100644
index 0000000..0f6ab50
--- /dev/null
+++ b/tests/reference_data/spatial_gs_energy.json
@@ -0,0 +1,38 @@
+{
+  "mp": {
+    "0": {
+      "r": "2 {f^{i_{\\alpha}}_{i_{\\alpha}}}",
+      "u": "{f^{i_{\\alpha}}_{i_{\\alpha}}} + {f^{i_{\\beta}}_{i_{\\beta}}}"
+    },
+    "1": {
+      "r": "- 2 {v^{i_{\\alpha}i_{\\alpha}}_{j_{\\alpha}j_{\\alpha}}} + {v^{i_{\\alpha}j_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}}",
+      "u": "- {v^{i_{\\alpha}i_{\\alpha}}_{i_{\\beta}i_{\\beta}}} - \\frac{{v^{i_{\\alpha}i_{\\alpha}}_{j_{\\alpha}j_{\\alpha}}}}{2} + \\frac{{v^{i_{\\alpha}j_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}}}{2} - \\frac{{v^{i_{\\beta}i_{\\beta}}_{j_{\\beta}j_{\\beta}}}}{2} + \\frac{{v^{i_{\\beta}j_{\\beta}}_{i_{\\beta}j_{\\beta}}}}{2}"
+    },
+    "2": {
+      "r": "- \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}a_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{2} + \\frac{3 {t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}b_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{2}",
+      "u": "- {t1^{a_{\\alpha}a_{\\beta}}_{i_{\\alpha}i_{\\beta}}} {v^{i_{\\alpha}a_{\\alpha}}_{i_{\\beta}a_{\\beta}}} - \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}a_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{4} + \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}b_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{4} - \\frac{{t1^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {v^{i_{\\beta}a_{\\beta}}_{j_{\\beta}b_{\\beta}}}}{4} + \\frac{{t1^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {v^{i_{\\beta}b_{\\beta}}_{j_{\\beta}a_{\\beta}}}}{4}"
+    },
+    "3": {
+      "r": "- \\frac{3 {t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}a_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{2} + \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}b_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{2}",
+      "u": "- {t2^{a_{\\alpha}a_{\\beta}}_{i_{\\alpha}i_{\\beta}}} {v^{i_{\\alpha}a_{\\alpha}}_{i_{\\beta}a_{\\beta}}} - \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}a_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{4} + \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}b_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{4} - \\frac{{t2^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {v^{i_{\\beta}a_{\\beta}}_{j_{\\beta}b_{\\beta}}}}{4} + \\frac{{t2^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {v^{i_{\\beta}b_{\\beta}}_{j_{\\beta}a_{\\beta}}}}{4}"
+    }
+  },
+  "re": {
+    "0": {
+      "r": "2 {f^{i_{\\alpha}}_{i_{\\alpha}}} - 2 {v^{i_{\\alpha}i_{\\alpha}}_{j_{\\alpha}j_{\\alpha}}} + {v^{i_{\\alpha}j_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}}",
+      "u": "{f^{i_{\\alpha}}_{i_{\\alpha}}} + {f^{i_{\\beta}}_{i_{\\beta}}} - {v^{i_{\\alpha}i_{\\alpha}}_{i_{\\beta}i_{\\beta}}} - \\frac{{v^{i_{\\alpha}i_{\\alpha}}_{j_{\\alpha}j_{\\alpha}}}}{2} + \\frac{{v^{i_{\\alpha}j_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}}}{2} - \\frac{{v^{i_{\\beta}i_{\\beta}}_{j_{\\beta}j_{\\beta}}}}{2} + \\frac{{v^{i_{\\beta}j_{\\beta}}_{i_{\\beta}j_{\\beta}}}}{2}"
+    },
+    "1": {
+      "r": "0",
+      "u": "0"
+    },
+    "2": {
+      "r": "- \\frac{3 {t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}a_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{2} + \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}b_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{2}",
+      "u": "- {t1^{a_{\\alpha}a_{\\beta}}_{i_{\\alpha}i_{\\beta}}} {v^{i_{\\alpha}a_{\\alpha}}_{i_{\\beta}a_{\\beta}}} - \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}a_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{4} + \\frac{{t1^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}b_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{4} - \\frac{{t1^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {v^{i_{\\beta}a_{\\beta}}_{j_{\\beta}b_{\\beta}}}}{4} + \\frac{{t1^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {v^{i_{\\beta}b_{\\beta}}_{j_{\\beta}a_{\\beta}}}}{4}"
+    },
+    "3": {
+      "r": "2 {t2^{a_{\\alpha}}_{i_{\\alpha}}} {f^{i_{\\alpha}}_{a_{\\alpha}}} - \\frac{3 {t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}a_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{2} + \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}b_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{2}",
+      "u": "{t2^{a_{\\alpha}}_{i_{\\alpha}}} {f^{i_{\\alpha}}_{a_{\\alpha}}} + {t2^{a_{\\beta}}_{i_{\\beta}}} {f^{i_{\\beta}}_{a_{\\beta}}} - {t2^{a_{\\alpha}a_{\\beta}}_{i_{\\alpha}i_{\\beta}}} {v^{i_{\\alpha}a_{\\alpha}}_{i_{\\beta}a_{\\beta}}} - \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}a_{\\alpha}}_{j_{\\alpha}b_{\\alpha}}}}{4} + \\frac{{t2^{a_{\\alpha}b_{\\alpha}}_{i_{\\alpha}j_{\\alpha}}} {v^{i_{\\alpha}b_{\\alpha}}_{j_{\\alpha}a_{\\alpha}}}}{4} - \\frac{{t2^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {v^{i_{\\beta}a_{\\beta}}_{j_{\\beta}b_{\\beta}}}}{4} + \\frac{{t2^{a_{\\beta}b_{\\beta}}_{i_{\\beta}j_{\\beta}}} {v^{i_{\\beta}b_{\\beta}}_{j_{\\beta}a_{\\beta}}}}{4}"
+    }
+  }
+}
\ No newline at end of file
diff --git a/tests/resolution_of_identity_test.py b/tests/resolution_of_identity_test.py
new file mode 100644
index 0000000..fb745d4
--- /dev/null
+++ b/tests/resolution_of_identity_test.py
@@ -0,0 +1,32 @@
+from adcgen.spatial_orbitals import transform_to_spatial_orbitals
+from adcgen.resolution_of_identity import apply_resolution_of_identity
+from adcgen.simplify import simplify
+from adcgen.expression import ExprContainer
+
+from sympy import S
+
+import pytest
+
+
+class TestResolutionOfIdentity():
+
+    @pytest.mark.parametrize('variant', ['mp', 're'])
+    @pytest.mark.parametrize('order', [0, 1, 2, 3])
+    @pytest.mark.parametrize('restriction', ['r', 'u'])
+    @pytest.mark.parametrize('symmetry', ['sym', 'asym'])
+    def test_ri_gs_energy(self, variant, order, restriction, symmetry,
+                          cls_instances, reference_data):
+        # load the reference data
+        ref = reference_data['ri_gs_energy'][variant][order]
+        ref = ref[restriction][symmetry]
+        ref.make_real()
+        # transform restriction and symmetry to bool
+        restricted = restriction == 'r'
+        # compute the energy
+        e = cls_instances[variant]['gs'].energy(order)
+        expr = ExprContainer(e, real=True)
+
+        sp_expr = transform_to_spatial_orbitals(expr, '', '', restricted)
+        ri_expr = apply_resolution_of_identity(sp_expr, symmetry)
+
+        assert simplify(ri_expr - ref).substitute_contracted().inner is S.Zero
diff --git a/tests/spatial_orbitals_test.py b/tests/spatial_orbitals_test.py
index aa95aba..405f8e4 100644
--- a/tests/spatial_orbitals_test.py
+++ b/tests/spatial_orbitals_test.py
@@ -14,6 +14,8 @@
 from sympy import Add, S, Rational, sympify
 from sympy.physics.secondquant import F, Fd
 
+import pytest
+
 
 class TestExpandAntiSymEri:
     def test_no_eri(self):
@@ -338,3 +340,26 @@ def test_t1_2(self):
         target = restricted.provided_target_idx
         assert target is not None
         assert set(target) == {i, a}
+
+
+class TestSpatialGroundstateEnergy:
+
+    @pytest.mark.parametrize('variant', ['mp', 're'])
+    @pytest.mark.parametrize('order', [0, 1, 2])
+    @pytest.mark.parametrize('restriction', ['r', 'u'])
+    def test_spatial_gs_energy(self, variant, order, restriction,
+                               cls_instances, reference_data):
+        # load the reference data
+        ref = reference_data['spatial_gs_energy'][variant][order][restriction]
+        # transform restriction to bool
+        restricted = restriction == 'r'
+        # compute the energy
+        e = cls_instances[variant]['gs'].energy(order)
+        expr = ExprContainer(e, real=True)
+
+        sp_expr = transform_to_spatial_orbitals(expr, '', '', restricted)
+        ref.make_real()
+        ref.add_bra_ket_sym(braket_sym_tensors=tensor_names.coulomb)
+
+        assert (simplify(sp_expr - ref).substitute_contracted().inner
+                is S.Zero)