diff --git a/src/qforte/adapters/molecule_adapters.py b/src/qforte/adapters/molecule_adapters.py
index 98692417..17df086b 100644
--- a/src/qforte/adapters/molecule_adapters.py
+++ b/src/qforte/adapters/molecule_adapters.py
@@ -48,9 +48,6 @@ def create_psi_mol(**kwargs):
     #[2, 0, 0, 1] for an A1A1B2 determinant.
     kwargs.setdefault('scf_docc', None)
 
-    #Reference wavefunction from previous casscf calculation
-    kwargs.setdefault('ref_wfn', None)
-
     #Tuple containing restricted occupied, active, and restricted virtual irrep indices
     #e.g. H4 in D2h symmetry with B2u and B3u orbitals as the active space is
     #([1,0,...0],[0,...,1,1],[0,1,0,...0])
@@ -128,9 +125,9 @@ def create_psi_mol(**kwargs):
         psi4.set_options({'mcscf_r_convergence': 1e-12})  
         psi4.set_options({'mcscf_maxiter': 1000})  
         psi4.set_options({'mcscf_diis_start': 50})   
-        E_casscf, p4_wfn = psi4.energy('casscf', return_wfn=True, ref_wfn = kwargs['ref_wfn'])   
+        E_casscf, p4_wfn = psi4.energy('casscf', return_wfn=True)   
         p4_Escf = None
-
+        print(f"CASSCF Energy: {E_casscf}")
     else:  
         p4_Escf, p4_wfn = psi4.energy('SCF', return_wfn=True)
         
diff --git a/src/qforte/expansions/manifolds.py b/src/qforte/expansions/manifolds.py
index 0e0dd77e..40717269 100644
--- a/src/qforte/expansions/manifolds.py
+++ b/src/qforte/expansions/manifolds.py
@@ -86,11 +86,13 @@ def sa_single(ref, p, q, mult = 1):
     """
     Get a specific, spin-adapted single excitation a_p^q|ref>.  Mult determines the singlet or triplet CSF.
     """
+ 
+    pa, pb, qa, qb = 2*p, 2*p + 1, 2*q, 2*q + 1
     diff = [0]*len(ref)
-    diff[p] = -1
-    diff[q] = 1
+    diff[pa] = -1
+    diff[qa] = 1
     new_det = [ref[i] + diff[i] for i in range(len(ref))]
-    pa, pb, qa, qb = 2*p, 2*p + 1, 2*q, 2*q + 1
+    
     U = qf.Circuit()
     do_j = []
     if mult == 1:    
@@ -99,13 +101,14 @@ def sa_single(ref, p, q, mult = 1):
                 do_j.append(2*k) 
                 if 2*k+1 not in [pb, qb]:
                     do_j.append(2*k + 1) 
-            for j in do_j:
-                U.add(qf.gate('X', j, j)) 
-            U.add(qf.gate("H", pa, pa))
-            U.add(qf.gate("CNOT", pb, pa))
-            U.add(qf.gate("X", pa, pa))
-            U.add(qf.gate("CNOT", qa, pa))
-            U.add(qf.gate("CNOT", qb, pa))
+        for j in do_j:
+            U.add(qf.gate('X', j, j)) 
+        U.add(qf.gate("H", pa, pa))
+        U.add(qf.gate("CNOT", pb, pa))
+        U.add(qf.gate("X", pa, pa))
+        U.add(qf.gate("CNOT", qa, pa))
+        U.add(qf.gate("CNOT", qb, pa))
+            
     elif mult == 3:
         for k in range(int(len(ref)/2)):
             if new_det[2 * k] == 1 or new_det[2*k + 1] == 1: