diff --git a/pyerrors/input/hadrons.py b/pyerrors/input/hadrons.py
index 1f21c3e0..00b93453 100644
--- a/pyerrors/input/hadrons.py
+++ b/pyerrors/input/hadrons.py
@@ -63,7 +63,7 @@ def read_meson_hd5(path, filestem, ens_id, meson='meson_0', tree='meson', idl=No
         for outputs of the Meson module. Can be altered to read input
         from other modules with similar structures.
     idl : range
-        If specified only conifgurations in the given range are read in.
+        If specified only configurations in the given range are read in.
     """
 
     files, idx = _get_files(path, filestem, idl)
@@ -134,12 +134,15 @@ def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
     """Read hadrons ExternalLeg hdf5 file and output an array of CObs
 
     Parameters
-    -----------------
-    path -- path to the files to read
-    filestem -- namestem of the files to read
-    ens_id -- name of the ensemble, required for internal bookkeeping
+    ----------
+    path : str
+        path to the files to read
+    filestem : str
+        namestem of the files to read
+    ens_id : str
+        name of the ensemble, required for internal bookkeeping
     idl : range
-        If specified only conifgurations in the given range are read in.
+        If specified only configurations in the given range are read in.
     """
 
     files, idx = _get_files(path, filestem, idl)
@@ -171,12 +174,15 @@ def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
     """Read hadrons Bilinear hdf5 file and output an array of CObs
 
     Parameters
-    -----------------
-    path -- path to the files to read
-    filestem -- namestem of the files to read
-    ens_id -- name of the ensemble, required for internal bookkeeping
+    ----------
+    path : str
+        path to the files to read
+    filestem : str
+        namestem of the files to read
+    ens_id : str
+        name of the ensemble, required for internal bookkeeping
     idl : range
-        If specified only conifgurations in the given range are read in.
+        If specified only configurations in the given range are read in.
     """
 
     files, idx = _get_files(path, filestem, idl)
@@ -216,3 +222,62 @@ def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
         result_dict[key] = Npr_matrix(matrix.swapaxes(1, 2).reshape((12, 12), order='F'), mom_in=mom_in, mom_out=mom_out)
 
     return result_dict
+
+
+def read_Fourquark_hd5(path, filestem, ens_id, idl=None):
+    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
+
+    Parameters
+    ----------
+    path : str
+        path to the files to read
+    filestem : str
+        namestem of the files to read
+    ens_id : str
+        name of the ensemble, required for internal bookkeeping
+    idl : range
+        If specified only configurations in the given range are read in.
+    """
+
+    files, idx = _get_files(path, filestem, idl)
+
+    mom_in = None
+    mom_out = None
+
+    corr_data = {}
+
+    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
+
+    for hd5_file in files:
+        file = h5py.File(path + '/' + hd5_file, "r")
+
+        for i in range(1):
+            name = file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8') + '_' + file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8')
+            if name not in corr_data:
+                corr_data[name] = []
+            raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
+            corr_data[name].append(raw_data)
+            if mom_in is None:
+                mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(' '), dtype=int)
+            if mom_out is None:
+                mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(' '), dtype=int)
+
+        file.close()
+
+    result_dict = {}
+
+    for key, data in corr_data.items():
+        local_data = np.array(data)
+
+        rolled_array = np.moveaxis(local_data, 0, 8)
+
+        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
+        for index in np.ndindex(rolled_array.shape[:-1]):
+            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
+            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
+            matrix[index] = CObs(real, imag)
+
+        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
+        # result_dict[key] = Npr_matrix(matrix.swapaxes(1, 2).reshape((12, 12), order='F'), mom_in=mom_in, mom_out=mom_out)
+
+    return result_dict
diff --git a/pyerrors/linalg.py b/pyerrors/linalg.py
index fb121e6a..dbbde944 100644
--- a/pyerrors/linalg.py
+++ b/pyerrors/linalg.py
@@ -186,49 +186,49 @@ def jack_matmul(*operands):
     For large matrices this is considerably faster compared to matmul.
     """
 
-    if any(isinstance(o[0, 0], CObs) for o in operands):
-        name = operands[0][0, 0].real.names[0]
-        idl = operands[0][0, 0].real.idl[name]
+    if any(isinstance(o.flat[0], CObs) for o in operands):
+        name = operands[0].flat[0].real.names[0]
+        idl = operands[0].flat[0].real.idl[name]
 
         def _exp_to_jack(matrix):
             base_matrix = np.empty_like(matrix)
-            for (n, m), entry in np.ndenumerate(matrix):
-                base_matrix[n, m] = entry.real.export_jackknife() + 1j * entry.imag.export_jackknife()
+            for index, entry in np.ndenumerate(matrix):
+                base_matrix[index] = entry.real.export_jackknife() + 1j * entry.imag.export_jackknife()
             return base_matrix
 
         def _imp_from_jack(matrix):
             base_matrix = np.empty_like(matrix)
-            for (n, m), entry in np.ndenumerate(matrix):
-                base_matrix[n, m] = CObs(import_jackknife(entry.real, name, [idl]),
-                                         import_jackknife(entry.imag, name, [idl]))
+            for index, entry in np.ndenumerate(matrix):
+                base_matrix[index] = CObs(import_jackknife(entry.real, name, [idl]),
+                                          import_jackknife(entry.imag, name, [idl]))
             return base_matrix
 
         r = _exp_to_jack(operands[0])
         for op in operands[1:]:
-            if isinstance(op[0, 0], CObs):
+            if isinstance(op.flat[0], CObs):
                 r = r @ _exp_to_jack(op)
             else:
                 r = r @ op
         return _imp_from_jack(r)
     else:
-        name = operands[0][0, 0].names[0]
-        idl = operands[0][0, 0].idl[name]
+        name = operands[0].flat[0].names[0]
+        idl = operands[0].flat[0].idl[name]
 
         def _exp_to_jack(matrix):
             base_matrix = np.empty_like(matrix)
-            for (n, m), entry in np.ndenumerate(matrix):
-                base_matrix[n, m] = entry.export_jackknife()
+            for index, entry in np.ndenumerate(matrix):
+                base_matrix[index] = entry.export_jackknife()
             return base_matrix
 
         def _imp_from_jack(matrix):
             base_matrix = np.empty_like(matrix)
-            for (n, m), entry in np.ndenumerate(matrix):
-                base_matrix[n, m] = import_jackknife(entry, name, [idl])
+            for index, entry in np.ndenumerate(matrix):
+                base_matrix[index] = import_jackknife(entry, name, [idl])
             return base_matrix
 
         r = _exp_to_jack(operands[0])
         for op in operands[1:]:
-            if isinstance(op[0, 0], Obs):
+            if isinstance(op.flat[0], Obs):
                 r = r @ _exp_to_jack(op)
             else:
                 r = r @ op