diff --git a/.DS_Store b/.DS_Store
index 2ba9821..97f546b 100644
Binary files a/.DS_Store and b/.DS_Store differ
diff --git a/symmer/.DS_Store b/symmer/.DS_Store
index 1b6e71e..b80ade3 100644
Binary files a/symmer/.DS_Store and b/symmer/.DS_Store differ
diff --git a/symmer/__init__.py b/symmer/__init__.py
index 5252ee9..4f59688 100644
--- a/symmer/__init__.py
+++ b/symmer/__init__.py
@@ -1,4 +1,6 @@
 """Main init for package."""
+import warnings
+warnings.filterwarnings('ignore', module='cotengra')
 from symmer.process_handler import process
 from symmer.operators  import PauliwordOp, QuantumState
 from symmer.projection import QubitTapering, ContextualSubspace, QubitSubspaceManager
\ No newline at end of file
diff --git a/symmer/operators/base.py b/symmer/operators/base.py
index 5eaaced..c1fb59e 100644
--- a/symmer/operators/base.py
+++ b/symmer/operators/base.py
@@ -21,6 +21,11 @@
 from qiskit.quantum_info import SparsePauliOp
 warnings.simplefilter('always', UserWarning)
 
+from qiskit._accelerate.sparse_pauli_op import (
+    ZXPaulis,
+    to_matrix_sparse,
+)
+
 from numba.core.errors import NumbaDeprecationWarning, NumbaPendingDeprecationWarning
 warnings.simplefilter('ignore', category=NumbaDeprecationWarning)
 warnings.simplefilter('ignore', category=NumbaPendingDeprecationWarning)
@@ -1464,33 +1469,45 @@ def to_sparse_matrix(self) -> csr_matrix:
         if self.n_qubits == 0:
             return csr_matrix(self.coeff_vec)
 
-        if self.n_qubits>15:
-            from symmer.utils import get_sparse_matrix_large_pauliwordop
-            sparse_matrix = get_sparse_matrix_large_pauliwordop(self)
-            return sparse_matrix
-        else:
-            x_int = binary_array_to_int(self.X_block).reshape(-1, 1)
-            z_int = binary_array_to_int(self.Z_block).reshape(-1, 1)
+        # if self.n_qubits>15:
+        #     from symmer.utils import get_sparse_matrix_large_pauliwordop
+        #     sparse_matrix = get_sparse_matrix_large_pauliwordop(self)
+        #     return sparse_matrix
+        # else:
+        #     x_int = binary_array_to_int(self.X_block).reshape(-1, 1)
+        #     z_int = binary_array_to_int(self.Z_block).reshape(-1, 1)
 
-            Y_number = np.sum(np.bitwise_and(self.X_block, self.Z_block), axis=1)
-            global_phase = (-1j) ** Y_number
+        #     Y_number = np.sum(np.bitwise_and(self.X_block, self.Z_block), axis=1)
+        #     global_phase = (-1j) ** Y_number
 
-            dimension = 2 ** self.n_qubits
-            row_ind = np.repeat(np.arange(dimension).reshape(1, -1), self.X_block.shape[0], axis=0)
-            col_ind = np.bitwise_xor(row_ind, x_int)
+        #     dimension = 2 ** self.n_qubits
+        #     row_ind = np.repeat(np.arange(dimension).reshape(1, -1), self.X_block.shape[0], axis=0)
+        #     col_ind = np.bitwise_xor(row_ind, x_int)
 
-            row_inds_and_Zint = np.bitwise_and(row_ind, z_int)
-            vals = global_phase.reshape(-1, 1) * (-1) ** (
-                        count1_in_int_bitstring(row_inds_and_Zint) % 2)  # .astype(complex))
+        #     row_inds_and_Zint = np.bitwise_and(row_ind, z_int)
+        #     vals = global_phase.reshape(-1, 1) * (-1) ** (
+        #                 count1_in_int_bitstring(row_inds_and_Zint) % 2)  # .astype(complex))
 
-            values_and_coeff = np.einsum('ij,i->ij', vals, self.coeff_vec)
+        #     values_and_coeff = np.einsum('ij,i->ij', vals, self.coeff_vec)
 
-            sparse_matrix = csr_matrix(
-                (values_and_coeff.flatten(), (row_ind.flatten(), col_ind.flatten())),
-                shape=(dimension, dimension),
-                dtype=complex
-            )
-            return sparse_matrix
+        #     sparse_matrix = csr_matrix(
+        #         (values_and_coeff.flatten(), (row_ind.flatten(), col_ind.flatten())),
+        #         shape=(dimension, dimension),
+        #         dtype=complex
+        #     )
+        #     return sparse_matrix
+        
+        phase = np.zeros(self.n_terms, dtype=np.uint8)
+        zx = ZXPaulis(
+                    self.X_block[:,::-1],
+                    self.Z_block[:,::-1],
+                    phase,
+                    self.coeff_vec,
+                )
+        
+        data, indices, indptr = to_matrix_sparse(zx, force_serial=False)
+        side = 1 << self.n_qubits
+        return csr_matrix((data, indices, indptr), shape=(side, side))
 
     def conjugate_op(self, R: 'PauliwordOp') -> 'PauliwordOp':
         """
diff --git a/symmer/process_handler.py b/symmer/process_handler.py
index aeccc93..f438eba 100644
--- a/symmer/process_handler.py
+++ b/symmer/process_handler.py
@@ -5,14 +5,18 @@
 from ray import remote, put, get
 from multiprocessing import Process, Queue, set_start_method
 
-if sys.platform.lower() in ['linux', 'darwin']:
+if sys.platform.lower() in ['linux', 'darwin', 'linux2']:
     set_start_method('fork', force = True)
 else:
     set_start_method('spawn', force = True)
 
 class ProcessHandler:
 
-    method  = 'ray'
+    if sys.platform.lower() in ['linux', 'darwin', 'linux2']:
+        method  = 'mp'
+    else:
+        method  = 'ray'
+    
     verbose = False
 
     def __init__(self):