diff --git a/adcgen/__init__.py b/adcgen/__init__.py
index 4fc56f7..cc81810 100644
--- a/adcgen/__init__.py
+++ b/adcgen/__init__.py
@@ -3,7 +3,7 @@
 from .eri_orbenergy import EriOrbenergy
 from .expression import ExprContainer, import_from_sympy_latex
 from .factor_intermediates import factor_intermediates
-from .func import evaluate_deltas, wicks
+from .func import evaluate_deltas
 from .generate_code import (generate_code, optimize_contractions, Contraction,
                             unoptimized_contraction)
 from .groundstate import GroundState
@@ -21,6 +21,7 @@
                             NonSymmetricTensor, KroneckerDelta, SymbolicTensor)
 from .tensor_names import tensor_names
 from .resolution_of_identity import apply_resolution_of_identity
+from .wicks import wicks
 from . import sort_expr as sort
 
 
diff --git a/adcgen/func.py b/adcgen/func.py
index 6f26c88..71b23c8 100644
--- a/adcgen/func.py
+++ b/adcgen/func.py
@@ -1,15 +1,10 @@
 from collections.abc import Sequence
 import itertools
 
-from sympy.physics.secondquant import (
-    F, Fd, FermionicOperator, NO
-)
-from sympy import S, Add, Expr, Mul
+from sympy import Add, Expr, Mul
 
-from .expression import ExprContainer
 from .misc import Inputerror
-from .rules import Rules
-from .indices import Index, Indices, get_symbols
+from .indices import Index, get_symbols
 from .sympy_objects import KroneckerDelta
 
 
@@ -127,173 +122,3 @@ def evaluate_deltas(
         return expr
     else:
         return expr
-
-
-def wicks(expr: Expr, rules: Rules | None = None,
-          simplify_kronecker_deltas: bool = False) -> Expr:
-    """
-    Evaluates Wicks theorem in the provided expression only returning fully
-    contracted contributions.
-    Adapted from the implementation in 'sympy.physics.secondquant'.
-
-    Parameters
-    ----------
-    expr: Expr
-        Expression containing the second quantized operator strings to
-        evaluate. This function expects plain sympy objects (Add/Mul/...)
-        and no container class.
-    rules : Rules, optional
-        Rules that are applied to the result before returning, e.g., in the
-        context of RE not all tensor blocks might be allowed in the result.
-    simplify_kronecker_deltas : bool, optional
-        If set, the KroneckerDeltas generated through the contractions
-        will be evaluated before returning.
-    """
-    assert isinstance(expr, Expr)
-    # normal ordered operator string has to evaluate to zero
-    # and a single second quantized operator can not be contracted
-    if isinstance(expr, (NO, FermionicOperator)):
-        return S.Zero
-
-    # break up any NO-objects, and evaluate commutators
-    expr = expr.doit(wicks=True).expand()
-    assert isinstance(expr, Expr)
-
-    if isinstance(expr, Add):
-        return Add(*(
-            wicks(term, rules=rules,
-                  simplify_kronecker_deltas=simplify_kronecker_deltas)
-            for term in expr.args
-        ))
-    elif not isinstance(expr, Mul):
-        # nether Add, Mul, NO, F, Fd -> maybe a number or tensor
-        return expr
-    # -> we have a Mul object
-    # we don't want to mess around with commuting part of Mul
-    # so we factorize it out before starting recursion
-    c_part: list[Expr] = []
-    op_string: list[FermionicOperator] = []
-    for factor in expr.args:
-        if factor.is_commutative:
-            c_part.append(factor)
-        else:
-            assert isinstance(factor, FermionicOperator)
-            op_string.append(factor)
-
-    if (n := len(op_string)) == 0:  # no operators
-        result = expr
-    elif n == 1:  # a single operator
-        return S.Zero
-    else:  # at least 2 operators
-        result = _contract_operator_string(op_string)
-        result = Mul(*c_part, result).expand()
-        assert isinstance(result, Expr)
-        if simplify_kronecker_deltas:
-            result = evaluate_deltas(result)
-
-    # apply rules to the result
-    if rules is None:
-        return result
-    assert isinstance(rules, Rules)
-    return rules.apply(ExprContainer(result)).inner
-
-
-def _contract_operator_string(op_string: list[FermionicOperator]) -> Expr:
-    """
-    Contracts the operator string only returning fully contracted
-    contritbutions.
-    Adapted from 'sympy.physics.secondquant'.
-    """
-    # check that we can get a fully contracted contribution
-    if not _has_fully_contracted_contribution(op_string):
-        return S.Zero
-
-    result = []
-    for i in range(1, len(op_string)):
-        c = _contraction(op_string[0], op_string[i])
-        if c is S.Zero:
-            continue
-        if not i % 2:  # introduce -1 for swapping operators
-            c *= S.NegativeOne
-
-        if len(op_string) - 2 > 0:  # at least one operator left
-            # remove the contracted operators from the string and recurse
-            remaining = op_string[1:i] + op_string[i+1:]
-            result.append(Mul(c, _contract_operator_string(remaining)))
-        else:  # no operators left
-            result.append(c)
-    return Add(*result)
-
-
-def _contraction(p: FermionicOperator, q: FermionicOperator) -> Expr:
-    """
-    Evaluates the contraction between two sqcond quantized fermionic
-    operators.
-    Adapted from 'sympy.physics.secondquant'.
-    """
-    assert isinstance(p, FermionicOperator)
-    assert isinstance(q, FermionicOperator)
-
-    p_idx, q_idx = p.args[0], q.args[0]
-    assert isinstance(p_idx, Index) and isinstance(q_idx, Index)
-    if p_idx.spin or q_idx.spin:
-        raise NotImplementedError("Contraction not implemented for indices "
-                                  "with spin.")
-    # get the space and ensure we have no unexpected space
-    space_p, space_q = p_idx.space[0], q_idx.space[0]
-    assert space_p in ["o", "v", "g"] and space_q in ["o", "v", "g"]
-
-    if isinstance(p, F) and isinstance(q, Fd):
-        if space_p == "o" or space_q == "o":
-            res = S.Zero
-        elif space_p == "v" or space_q == "v":
-            res = KroneckerDelta(p_idx, q_idx)
-        else:
-            res = Mul(
-                KroneckerDelta(p_idx, q_idx),
-                KroneckerDelta(q_idx, Index('a', above_fermi=True))
-            )
-    elif isinstance(p, Fd) and isinstance(q, F):
-        if space_p == "v" or space_q == "v":
-            res = S.Zero
-        elif space_p == "o" or space_q == "o":
-            res = KroneckerDelta(p_idx, q_idx)
-        else:
-            res = Mul(
-                KroneckerDelta(p_idx, q_idx),
-                KroneckerDelta(q_idx, Index('i', below_fermi=True))
-            )
-    else:  # vanish if 2xAnnihilator or 2xCreator
-        res = S.Zero
-    assert isinstance(res, Expr)
-    return res
-
-
-def _has_fully_contracted_contribution(op_string: list[FermionicOperator]
-                                       ) -> bool:
-    """
-    Takes a list of second quantized operators and checks whether a
-    non-vanishing fully contracted contribution can exist.
-    """
-    if len(op_string) % 2:  # odd number of operators
-        return False
-    # count the number of creation and annihilation operators per space
-    create = {space: 0 for space in Indices.base}
-    annihilate = {space: 0 for space in Indices.base}
-    for op in op_string:
-        if isinstance(op, Fd):
-            counter = create
-        else:
-            counter = annihilate
-        idx = op.args[0]
-        assert isinstance(idx, Index)
-        counter[idx.space] += 1
-    # check that we have a matching amount of creation and annihilation
-    # operators
-    for space, n_create in create.items():
-        if space == "general":
-            continue
-        n_annihilate = annihilate[space] + annihilate["general"]
-        if n_create - n_annihilate > 0:
-            return False
-    return True
diff --git a/adcgen/groundstate.py b/adcgen/groundstate.py
index 4ee5128..6dbd2dd 100644
--- a/adcgen/groundstate.py
+++ b/adcgen/groundstate.py
@@ -4,7 +4,7 @@
 from sympy import Expr, Mul, Rational, S, latex
 
 from .expression import ExprContainer
-from .func import gen_term_orders, wicks
+from .func import gen_term_orders
 from .indices import Indices, n_ov_from_space
 from .logger import logger
 from .misc import cached_member, Inputerror, validate_input
@@ -12,6 +12,7 @@
 from .simplify import simplify
 from .sympy_objects import Amplitude
 from .tensor_names import tensor_names
+from .wicks import wicks
 
 
 class GroundState:
diff --git a/adcgen/intermediate_states.py b/adcgen/intermediate_states.py
index 2580ac6..2cc0e8f 100644
--- a/adcgen/intermediate_states.py
+++ b/adcgen/intermediate_states.py
@@ -5,7 +5,7 @@
 from sympy import Expr, Mul, Rational, S, latex, nsimplify, diff, symbols
 
 from .expression import ExprContainer
-from .func import gen_term_orders, wicks, evaluate_deltas
+from .func import gen_term_orders, evaluate_deltas
 from .groundstate import GroundState
 from .indices import (
     n_ov_from_space, repeated_indices, Indices, generic_indices_from_space
@@ -15,6 +15,7 @@
 from .simplify import simplify
 from .sympy_objects import Amplitude
 from .tensor_names import tensor_names
+from .wicks import wicks
 
 
 class IntermediateStates:
diff --git a/adcgen/properties.py b/adcgen/properties.py
index 41ed670..f395318 100644
--- a/adcgen/properties.py
+++ b/adcgen/properties.py
@@ -4,13 +4,14 @@
 from sympy import Add, Expr, S, sqrt, sympify
 
 from .expression import ExprContainer
-from .func import gen_term_orders, wicks
+from .func import gen_term_orders
 from .indices import n_ov_from_space, generic_indices_from_space
 from .intermediate_states import IntermediateStates
 from .misc import Inputerror, cached_member, transform_to_tuple, validate_input
 from .rules import Rules
 from .secular_matrix import SecularMatrix
 from .simplify import simplify
+from .wicks import wicks
 
 
 class Properties:
diff --git a/adcgen/secular_matrix.py b/adcgen/secular_matrix.py
index d374db1..9fdf2ed 100644
--- a/adcgen/secular_matrix.py
+++ b/adcgen/secular_matrix.py
@@ -4,7 +4,7 @@
 from sympy import Add, Expr, Mul, S, sqrt
 
 from .expression import ExprContainer
-from .func import gen_term_orders, wicks, evaluate_deltas
+from .func import gen_term_orders, evaluate_deltas
 from .groundstate import GroundState
 from .indices import (
     repeated_indices, Indices, generic_indices_from_space, n_ov_from_space
@@ -14,6 +14,7 @@
 from .operators import Operators
 from .rules import Rules
 from .simplify import simplify
+from .wicks import wicks
 
 
 class SecularMatrix:
diff --git a/adcgen/wicks.py b/adcgen/wicks.py
new file mode 100644
index 0000000..b1f415a
--- /dev/null
+++ b/adcgen/wicks.py
@@ -0,0 +1,276 @@
+from collections.abc import Sequence
+import itertools
+import math
+
+from sympy.physics.secondquant import F, Fd, FermionicOperator, NO
+from sympy import Add, Expr, Mul, S
+
+from .expression import ExprContainer
+from .func import evaluate_deltas
+from .indices import Index, Indices
+from .rules import Rules
+from .sympy_objects import KroneckerDelta
+
+
+def wicks(expr: Expr, rules: Rules | None = None,
+          simplify_kronecker_deltas: bool = False) -> Expr:
+    """
+    Evaluates Wicks theorem in the provided expression only returning fully
+    contracted contributions.
+    Adapted from the implementation in 'sympy.physics.secondquant'.
+
+    Parameters
+    ----------
+    expr: Expr
+        Expression containing the second quantized operator strings to
+        evaluate. This function expects plain sympy objects (Add/Mul/...)
+        and no container class.
+    rules : Rules, optional
+        Rules that are applied to the result before returning, e.g., in the
+        context of RE not all tensor blocks might be allowed in the result.
+    simplify_kronecker_deltas : bool, optional
+        If set, the :py:class:`KroneckerDelta` set generated through the
+        contractions will be evaluated before returning.
+    """
+    assert isinstance(expr, Expr)
+    # normal ordered operator string has to evaluate to zero
+    # and a single second quantized operator can not be contracted
+    if isinstance(expr, (NO, FermionicOperator)):
+        return S.Zero
+
+    # break up any NO-objects, and evaluate commutators
+    expr = expr.doit(wicks=True).expand()
+
+    if isinstance(expr, Add):
+        return Add(*(
+            wicks(term, rules=rules,
+                  simplify_kronecker_deltas=simplify_kronecker_deltas)
+            for term in expr.args
+        ))
+    elif isinstance(expr, Mul):
+        # we don't want to mess around with commuting part of Mul
+        # so we factorize it out before starting recursion
+        c_part: list[Expr] = []
+        op_string: list[FermionicOperator] = []
+        for factor in expr.args:
+            if factor.is_commutative:
+                c_part.append(factor)
+            else:
+                assert isinstance(factor, FermionicOperator)
+                op_string.append(factor)
+
+        if (n := len(op_string)) == 0:  # no operators
+            result = expr
+        elif n == 1:  # a single operator
+            return S.Zero
+        else:  # at least 2 operators
+            result = _contract_operator_string(
+                op_string=tuple(enumerate(op_string))
+            )
+            result = Mul(*c_part, result).expand()
+            if simplify_kronecker_deltas:
+                result = evaluate_deltas(result)
+    else:  # neither Add, Mul, NO, F, Fd -> maybe a number or tensor
+        result = expr
+
+    # apply rules to the result
+    if rules is not None:
+        assert isinstance(rules, Rules)
+        result = rules.apply(ExprContainer(result)).inner
+    return result
+
+
+def _contract_operator_string(
+        op_string: Sequence[tuple[int, FermionicOperator]],
+        n_total_operators: int | None = None,
+        contractions: Sequence[Expr] | None = None) -> Expr:
+    """
+    Contracts the operator string only returning fully contracted
+    contritbutions.
+
+    Parameters
+    ----------
+    op_string: Sequence[tuple[int, FermionicOperator]]
+        The list/tuple of second quantized operators to contract
+        along with their positions in the operator string
+    n_total_operators: int | None, optional
+        The total number of second quantized operators in the
+        operator string. By default, it is calculated from the
+        length of the contraction cache.
+    contractions: Sequence[Expr] | None, optional
+        Precomputed flattened array of contractions of operator
+        pairs. Only the upper triangular part
+        (excluding diagonal elements) of the pair matrix
+        is expected to be present in the cache: for N operators
+        a cache with N(N-1)/2 elements is expected. If not given
+        it will be calculated on the fly for the current
+        operator string.
+    """
+    # This function implements recursive depth first tree traversal
+
+    # check that we can get a fully contracted contribution
+    if not _has_fully_contracted_contribution(op_string):
+        return S.Zero
+    # we can precompute all relevant contractions as they can
+    # be reused often in an depth first graph traversal
+    # It is sufficient to precompute the upper triangle
+    #   1  2  3
+    # 1    x  x
+    # 2       x
+    # 3
+    # the flattened index of the element (i, j) that is part of
+    # the upper triangle (excluding diagonal elements) of a n x n matrix
+    # can be computed accoding to
+    # idx = (2*i*n - i**2 + 2*j - 3*i - 2) / 2
+    # using row index i, column index j and the number of elements
+    # along each dimension n.
+    # for a derivation see:
+    # math.stackexchange.com/questions/646117/how-to-find-a-function-mapping-matrix-indices
+    if contractions is None:
+        contractions = tuple(
+            _contraction(op1, op2) for (_, op1), (_, op2) in
+            itertools.combinations(op_string, 2)
+        )
+    # calculate the number of operators from the length of the
+    # contraction cache: n(n-1)/2 elements are in the cache
+    # required for the calculation of the flattened cache index
+    if n_total_operators is None:
+        n = 0.5 + math.sqrt(0.25 + 2 * len(contractions))
+        # if n is not integer, the number of entries in the
+        # contraction cache is not valid
+        assert n.is_integer()
+        n_total_operators = int(n)
+        del n
+    # left_pos and right_pos denote the positions in the
+    # original full operator string (before contracting any
+    # operators). This is required for the cache lookup.
+    # The loop index i denotes the position in the
+    # current version of the operator string potentially
+    # with previous contraction of other operators. Required to
+    # to determine the sign of a contraction
+    result = []
+    left_pos, _ = op_string[0]
+    for i in range(1, len(op_string)):
+        right_pos, _ = op_string[i]
+        # compute the flattened index for the (left_pos, right_pos)
+        # element of the upper triangular matrix as explained above
+        flattened_idx = (
+            2 * left_pos * n_total_operators - left_pos * left_pos
+            + 2 * right_pos - 3 * left_pos - 2
+        ) // 2
+        c = contractions[flattened_idx]
+        if c is S.Zero:
+            continue
+        if not i % 2:  # introduce -1 for swapping operators
+            c *= S.NegativeOne
+
+        if len(op_string) - 2 > 0:  # at least one operator left
+            # remove the contracted operators from the string and recurse
+            remaining = tuple(
+                ele for j, ele in
+                enumerate(op_string) if j != 0 and j != i
+            )
+            c *= _contract_operator_string(
+                op_string=remaining, n_total_operators=n_total_operators,
+                contractions=contractions
+            )
+            result.append(c)
+        else:  # no operators left
+            result.append(c)
+    return Add(*result)
+
+
+def _contraction(p: FermionicOperator, q: FermionicOperator) -> Expr:
+    """
+    Evaluates the contraction between two sqcond quantized fermionic
+    operators.
+    Adapted from 'sympy.physics.secondquant'.
+    """
+    assert isinstance(p, FermionicOperator)
+    assert isinstance(q, FermionicOperator)
+
+    p_idx, q_idx = p.args[0], q.args[0]
+    assert isinstance(p_idx, Index) and isinstance(q_idx, Index)
+    # in principle this also works for indices with spin since the
+    # KroneckerDelta handles spin correctly. However, mixed deltas
+    # \delta_{i j_{\alpha}} might not be evaluated by evaluate_deltas,
+    # because the preferred index would be i_{\alpha} - a new index
+    # that is not present in the expression.
+    # To avoid unevaluated deltas this check was introduced
+    if p_idx.spin or q_idx.spin:
+        raise NotImplementedError("Contraction not implemented for indices "
+                                  "with spin.")
+    # get the space and ensure we have no unexpected space
+    space_p, space_q = p_idx.space[0], q_idx.space[0]
+    assert space_p in ["o", "v", "g"] and space_q in ["o", "v", "g"]
+
+    if isinstance(p, F) and isinstance(q, Fd):
+        if space_p == "o" or space_q == "o":
+            res = S.Zero
+        elif space_p == "v" or space_q == "v":
+            res = KroneckerDelta(p_idx, q_idx)
+        else:
+            res = (
+                KroneckerDelta(p_idx, q_idx) *
+                KroneckerDelta(q_idx, Index('a', above_fermi=True))
+            )
+    elif isinstance(p, Fd) and isinstance(q, F):
+        if space_p == "v" or space_q == "v":
+            res = S.Zero
+        elif space_p == "o" or space_q == "o":
+            res = KroneckerDelta(p_idx, q_idx)
+        else:
+            res = (
+                KroneckerDelta(p_idx, q_idx) *
+                KroneckerDelta(q_idx, Index('i', below_fermi=True))
+            )
+    else:  # vanish if 2xAnnihilator or 2xCreator
+        res = S.Zero
+    return res
+
+
+def _has_fully_contracted_contribution(
+        op_string: Sequence[tuple[int, FermionicOperator]]) -> bool:
+    """
+    Takes a list of second quantized operators and their respective positions
+    in the operator string and checks whether a non-vanishing fully contracted
+    contribution can exist.
+    """
+    if len(op_string) % 2:  # odd number of operators
+        return False
+    # count the number of creation and annihilation operators per space
+    idx_spaces: tuple[str, ...] = tuple(Indices.base.keys())
+    creation: list[int] = [0 for _ in idx_spaces]
+    annihilation: list[int] = [0 for _ in idx_spaces]
+    for _, op in op_string:
+        if isinstance(op, Fd):
+            counter = creation
+        else:
+            counter = annihilation
+        idx = op.args[0]
+        assert isinstance(idx, Index)
+        counter[idx_spaces.index(idx.space)] += 1
+    # check that we have a matching amount of creation and annihilation
+    # operators
+    if sum(creation) != sum(annihilation):
+        return False
+    # remove the general operators from the counters
+    n_general_creation = creation.pop(idx_spaces.index("general"))
+    n_general_annihilation = annihilation.pop(idx_spaces.index("general"))
+    # ensure that we have a matching amount of creation and annihilation
+    # operators for each space accounting for general operators
+    for n_create, n_annihilate in zip(creation, annihilation, strict=True):
+        if n_create == n_annihilate:
+            continue
+        elif n_create > n_annihilate:
+            n_general_annihilation -= n_create - n_annihilate
+            if n_general_annihilation < 0:
+                return False
+        else:  # n_create < n_annihilate
+            n_general_creation -= n_annihilate - n_create
+            if n_general_creation < 0:
+                return False
+    # we have to have the same number of general creation and annihilation
+    # operators at this point. Otherwise sum(creation) != sum(annihilation)
+    # and we would have failed above
+    return True
diff --git a/tests/func_test.py b/tests/func_test.py
index 0c7bcdb..0747f16 100644
--- a/tests/func_test.py
+++ b/tests/func_test.py
@@ -1,12 +1,9 @@
-from adcgen.expression import ExprContainer
-from adcgen.func import (
-    _contraction, _contract_operator_string, wicks, evaluate_deltas
-)
-from adcgen.indices import Index, get_symbols
+from adcgen.func import evaluate_deltas
+from adcgen.indices import Index
 from adcgen.sympy_objects import KroneckerDelta
 
-from sympy import Add, Mul, S
-from sympy.physics.secondquant import F, Fd
+from sympy import Mul
+from sympy.physics.secondquant import F
 
 
 class TestEvaluateDeltas:
@@ -35,97 +32,3 @@ def test_with_target_idx(self):
         assert evaluate_deltas(test, j) == F(j)
         assert evaluate_deltas(test, p) == F(i)
         assert evaluate_deltas(test, (i, j)) == test
-
-
-class TestWicks:
-    def test_contraction(self):
-        i, j = Index("i", below_fermi=True), Index("j", below_fermi=True)
-        a, b = Index("a", above_fermi=True), Index("b", above_fermi=True)
-        p, q = Index("p"), Index("q")
-
-        assert _contraction(Fd(i), F(j)) == KroneckerDelta(i, j)
-        assert _contraction(F(i), Fd(j)) is S.Zero
-        assert _contraction(F(i), F(j)) is S.Zero
-        assert _contraction(Fd(i), Fd(j)) is S.Zero
-
-        assert _contraction(Fd(a), F(b)) is S.Zero
-        assert _contraction(F(a), Fd(b)) == KroneckerDelta(a, b)
-        assert _contraction(F(a), F(b)) is S.Zero
-        assert _contraction(Fd(a), Fd(b)) is S.Zero
-
-        zero = Add(
-            _contraction(Fd(p), F(q)),
-            -Mul(KroneckerDelta(p, q),
-                 KroneckerDelta(q, Index("i", below_fermi=True)))
-        )
-        zero = ExprContainer(zero, target_idx="").substitute_contracted()
-        assert zero.inner is S.Zero
-        zero = Add(
-            _contraction(F(p), Fd(q)),
-            -Mul(KroneckerDelta(p, q),
-                 KroneckerDelta(q, Index("a", above_fermi=True)))
-        )
-        zero = ExprContainer(zero, target_idx="").substitute_contracted()
-        assert zero.inner is S.Zero
-
-        assert _contraction(Fd(i), F(p)) == KroneckerDelta(i, p)
-        assert _contraction(F(p), Fd(i)) is S.Zero
-        assert _contraction(F(i), Fd(p)) is S.Zero
-        assert _contraction(Fd(p), F(i)) == KroneckerDelta(i, p)
-
-        assert _contraction(Fd(a), F(p)) is S.Zero
-        assert _contraction(F(p), Fd(a)) == KroneckerDelta(a, p)
-        assert _contraction(F(p), Fd(a)) == KroneckerDelta(a, p)
-        assert _contraction(Fd(a), F(p)) is S.Zero
-
-        assert _contraction(Fd(a), F(i)) is S.Zero
-        assert _contraction(F(i), Fd(a)) is S.Zero
-        assert _contraction(F(a), Fd(i)) is S.Zero
-        assert _contraction(Fd(i), F(a)) is S.Zero
-
-    def test_contract_operator_string(self):
-        i, j, a, b, p, q = get_symbols("ijabpq")
-
-        op_string = [Fd(i), F(a), Fd(b), F(j)]
-        ref = Mul(KroneckerDelta(i, j), KroneckerDelta(a, b))
-        assert _contract_operator_string(op_string) == ref
-
-        op_string = [Fd(i), F(a), Fd(p), F(q), Fd(b), F(j)]
-        ref = Add(
-            -Mul(KroneckerDelta(i, q), KroneckerDelta(a, b),
-                 KroneckerDelta(p, j)),
-            Mul(KroneckerDelta(i, j), KroneckerDelta(a, p),
-                KroneckerDelta(q, b)),
-            Mul(KroneckerDelta(i, j), KroneckerDelta(a, b),
-                KroneckerDelta(p, q),
-                KroneckerDelta(q, Index("i", below_fermi=True)))
-        )
-        res = _contract_operator_string(op_string).expand()
-        zero = ExprContainer(
-            Add(ref, -res), target_idx="ijab"
-        ).substitute_contracted()
-        assert zero.inner is S.Zero
-
-    def test_wicks(self):
-        i, j, a, b, p, q = get_symbols("ijabpq")
-
-        expr = Mul(
-            Fd(i), F(a), Fd(p), F(q), Fd(b), F(j), 2, KroneckerDelta(i, j)
-        )
-        ref = S.Zero
-        ref -= Mul(
-            KroneckerDelta(i, q), KroneckerDelta(a, b), KroneckerDelta(p, j)
-        )
-        ref += Mul(
-            KroneckerDelta(i, j), KroneckerDelta(a, p), KroneckerDelta(q, b)
-        )
-        ref += Mul(
-            KroneckerDelta(i, j), KroneckerDelta(a, b), KroneckerDelta(p, q),
-            KroneckerDelta(q, Index("i", below_fermi=True))
-        )
-        ref *= Mul(2, KroneckerDelta(i, j))
-        res = wicks(expr)
-        zero = ExprContainer(
-            ref.expand() - res, target_idx="ijab"
-        ).substitute_contracted()
-        assert zero.inner is S.Zero
diff --git a/tests/wicks_test.py b/tests/wicks_test.py
new file mode 100644
index 0000000..2b32e3a
--- /dev/null
+++ b/tests/wicks_test.py
@@ -0,0 +1,109 @@
+from adcgen.expression import ExprContainer
+from adcgen.indices import Index, get_symbols
+from adcgen.sympy_objects import KroneckerDelta
+from adcgen.wicks import _contraction, _contract_operator_string, wicks
+
+from sympy.physics.secondquant import F, Fd
+from sympy import S
+
+import pytest
+
+
+class TestWicks:
+    def test_contraction(self):
+        i, j = Index("i", below_fermi=True), Index("j", below_fermi=True)
+        a, b = Index("a", above_fermi=True), Index("b", above_fermi=True)
+        p, q = Index("p"), Index("q")
+
+        assert _contraction(Fd(i), F(j)) == KroneckerDelta(i, j)
+        assert _contraction(F(i), Fd(j)) is S.Zero
+        assert _contraction(F(i), F(j)) is S.Zero
+        assert _contraction(Fd(i), Fd(j)) is S.Zero
+
+        assert _contraction(Fd(a), F(b)) is S.Zero
+        assert _contraction(F(a), Fd(b)) == KroneckerDelta(a, b)
+        assert _contraction(F(a), F(b)) is S.Zero
+        assert _contraction(Fd(a), Fd(b)) is S.Zero
+
+        zero = (
+            _contraction(Fd(p), F(q))
+            - (KroneckerDelta(p, q) *
+               KroneckerDelta(q, Index("i", below_fermi=True)))
+        )
+        zero = ExprContainer(zero, target_idx="").substitute_contracted()
+        assert zero.inner is S.Zero
+        zero = (
+            _contraction(F(p), Fd(q))
+            - (KroneckerDelta(p, q) *
+               KroneckerDelta(q, Index("a", above_fermi=True)))
+        )
+        zero = ExprContainer(zero, target_idx="").substitute_contracted()
+        assert zero.inner is S.Zero
+
+        assert _contraction(Fd(i), F(p)) == KroneckerDelta(i, p)
+        assert _contraction(F(p), Fd(i)) is S.Zero
+        assert _contraction(F(i), Fd(p)) is S.Zero
+        assert _contraction(Fd(p), F(i)) == KroneckerDelta(i, p)
+
+        assert _contraction(Fd(a), F(p)) is S.Zero
+        assert _contraction(F(p), Fd(a)) == KroneckerDelta(a, p)
+        assert _contraction(F(p), Fd(a)) == KroneckerDelta(a, p)
+        assert _contraction(Fd(a), F(p)) is S.Zero
+
+        assert _contraction(Fd(a), F(i)) is S.Zero
+        assert _contraction(F(i), Fd(a)) is S.Zero
+        assert _contraction(F(a), Fd(i)) is S.Zero
+        assert _contraction(Fd(i), F(a)) is S.Zero
+
+    def test_contraction_spin(self):
+        i, j = get_symbols("ij", "ab")
+        with pytest.raises(NotImplementedError):
+            _contraction(Fd(i), F(j))
+
+    def test_contract_operator_string(self):
+        i, j, a, b, p, q = get_symbols("ijabpq")
+
+        op_string = [Fd(i), F(a), Fd(b), F(j)]
+        ref = KroneckerDelta(i, j) * KroneckerDelta(a, b)
+        assert _contract_operator_string(tuple(enumerate(op_string))) == ref
+
+        op_string = [Fd(i), F(a), Fd(p), F(q), Fd(b), F(j)]
+        ref = (
+            - (KroneckerDelta(i, q) * KroneckerDelta(a, b) *
+               KroneckerDelta(p, j))
+            + (KroneckerDelta(i, j) * KroneckerDelta(a, p) *
+               KroneckerDelta(q, b))
+            + (KroneckerDelta(i, j) * KroneckerDelta(a, b) *
+               KroneckerDelta(p, q) *
+               KroneckerDelta(q, Index("i", below_fermi=True)))
+        )
+        res = _contract_operator_string(tuple(enumerate(op_string))).expand()
+        zero = ExprContainer(
+            ref - res, target_idx="ijab"
+        ).substitute_contracted()
+        assert zero.inner is S.Zero
+
+    def test_wicks(self):
+        i, j, a, b, p, q = get_symbols("ijabpq")
+
+        expr = (
+            Fd(i) * F(a) * Fd(p) * F(q) * Fd(b) * F(j) * 2 *
+            KroneckerDelta(i, j)
+        )
+        ref = S.Zero
+        ref -= (
+            KroneckerDelta(i, q) * KroneckerDelta(a, b) * KroneckerDelta(p, j)
+        )
+        ref += (
+            KroneckerDelta(i, j) * KroneckerDelta(a, p) * KroneckerDelta(q, b)
+        )
+        ref += (
+            KroneckerDelta(i, j) * KroneckerDelta(a, b) * KroneckerDelta(p, q)
+            * KroneckerDelta(q, Index("i", below_fermi=True))
+        )
+        ref *= (2 * KroneckerDelta(i, j))
+        res = wicks(expr)
+        zero = ExprContainer(
+            ref.expand() - res, target_idx="ijab"
+        ).substitute_contracted()
+        assert zero.inner is S.Zero