diff --git a/qiskit_addon_dice_solver/dice_solver.py b/qiskit_addon_dice_solver/dice_solver.py
index 64f93a3..e561b40 100644
--- a/qiskit_addon_dice_solver/dice_solver.py
+++ b/qiskit_addon_dice_solver/dice_solver.py
@@ -302,22 +302,25 @@ def _call_dice(dice_dir: Path, mpirun_options: Sequence[str] | str | None) -> No
     dice_log_path = os.path.join(dice_dir, "dice_solver_logfile.log")
     if mpirun_options:
         if isinstance(mpirun_options, str):
-            mpirun_options = [mpirun_options]
+            mpirun_options = mpirun_options.split()
         dice_call = ["mpirun"] + list(mpirun_options) + [dice_path]
     else:
         dice_call = ["mpirun", dice_path]
 
     with open(dice_log_path, "w") as logfile:
-        try:
-            subprocess.run(
-                dice_call, cwd=dice_dir, stdout=logfile, stderr=logfile, check=True
-            )
-        except subprocess.CalledProcessError as e:
-            raise DiceExecutionError(
-                command=dice_call,
-                returncode=e.returncode,
-                log_path=dice_log_path,
-            ) from e
+        process = subprocess.run(
+            dice_call, cwd=dice_dir, stdout=logfile, stderr=logfile
+        )
+    rdm_path = dice_dir / "spin1RDM.0.0.txt"
+    # We check this manually because Dice is returning non-zero codes on successful executions,
+    # so we can't rely on the status code to tell us whether Dice executed succesfully. Unclear
+    # if this only happens on Dice/Riken branch.
+    if not os.path.exists(rdm_path):
+        raise DiceExecutionError(
+            command=dice_call,
+            returncode=process.returncode,
+            log_path=dice_log_path,
+        )
 
 
 def _write_input_files(
@@ -331,9 +334,10 @@ def _write_input_files(
     max_iter: int,
 ) -> None:
     """Prepare the Dice inputs in the specified directory."""
-    ### Move the FCI Dump to dice dir ###
-    shutil.copy(active_space_path, os.path.join(dice_dir, "fcidump.txt"))
-
+    ### Move the FCI Dump to dice dir if it's not already there. ###
+    dice_fci_path = Path(dice_dir) / "fcidump.txt"
+    if not os.path.exists(dice_fci_path):
+        shutil.copy(active_space_path, dice_fci_path)
     ### Write the input.dat ###
     num_elec = num_up + num_dn
     # Return only the lowest-energy state
diff --git a/test/test_n2.py b/test/test_n2.py
index 1439533..df5a549 100644
--- a/test/test_n2.py
+++ b/test/test_n2.py
@@ -15,14 +15,11 @@
 from pyscf import ao2mo, tools
 import numpy as np
 
-from sqd.utils.counts import generate_counts_uniform, counts_to_arrays
-from sqd.configuration_recovery import recover_configurations
-from qiskit_addon_dice_solver import solve_dice
-from sqd.subsampling import postselect_and_subsample
-from sqd.utils.fermion import (
-    bitstring_matrix_to_sorted_addresses,
-    flip_orbital_occupancies,
-)
+from qiskit_addon_sqd.counts import generate_counts_uniform, counts_to_arrays
+from qiskit_addon_sqd.configuration_recovery import recover_configurations
+from qiskit_addon_sqd.subsampling import postselect_and_subsample
+from qiskit_addon_sqd.fermion import flip_orbital_occupancies
+from qiskit_addon_dice_solver import solve_fermion
 
 
 # Specify molecule properties
@@ -77,35 +74,28 @@
             occupancies_bitwise,
             num_elec_a,
             num_elec_b,
-            # rand_seed=rand_seed,
         )
 
     # Throw out samples with incorrect hamming weight and create batches of subsamples.
     batches = postselect_and_subsample(
         bs_mat_tmp,
         probs_arr_tmp,
-        num_elec_a,
-        num_elec_b,
-        samples_per_batch,
-        n_batches,
-        # rand_seed=rand_seed,
+        hamming_right=num_elec_a,
+        hamming_left=num_elec_b,
+        samples_per_batch=samples_per_batch,
+        num_batches=n_batches,
     )
     # Run eigenstate solvers in a loop. This loop should be parallelized for larger problems.
     int_e = np.zeros(n_batches)
     int_occs = np.zeros((n_batches, 2 * num_orbitals))
-    for j in range(n_batches):
-        addresses = bitstring_matrix_to_sorted_addresses(
-            batches[j], open_shell=open_shell
-        )
-        energy_sci, wf_mags, avg_occs = solve_dice(
-            addresses,
-            active_space_path,
-            os.path.abspath(os.path.dirname(__file__)),
-            spin_sq=spin_sq,
-            max_davidson=max_davidson_cycles,
-            clean_working_dir=True,
-            mpirun_options=["-n", "8"],
+    for j, batch in enumerate(batches):
+        energy_sci, wf_mags, avg_occs = solve_fermion(
+            batch,
+            hcore,
+            eri,
+            mpirun_options=["-quiet", "-n", "8"],
         )
+        energy_sci += nuclear_repulsion_energy
         int_e[j] = energy_sci
         int_occs[j, :num_orbitals] = avg_occs[0]
         int_occs[j, num_orbitals:] = avg_occs[1]
@@ -119,4 +109,4 @@
     e_hist[i, :] = int_e
 
 print(f"Exact energy: -109.10288938")
-print(f"Estimated energy: {np.min(e_hist[-1]) + nuclear_repulsion_energy}")
+print(f"Estimated energy: {np.min(e_hist[-1])}")