diff --git a/src/iterDiag.jl b/src/iterDiag.jl
index f11ef44..2932006 100644
--- a/src/iterDiag.jl
+++ b/src/iterDiag.jl
@@ -115,21 +115,26 @@ function IterDiag(
     end
 
     operators = CreateProductOperator(create[1], operators, newSitesFlow[1])
+    totalNumOperator = sum([operators[("n", [site])] for site in currentSites])
 
     energyPerSite = Float64[]
     @showprogress for (step, hamlt) in enumerate(hamltFlow)
-        for (type, members, strength) in hamlt
+        @time for (type, members, strength) in hamlt
             hamltMatrix += strength * operators[(type, members)]
         end
+        println("Hamiltonian size = ", size(hamltMatrix))
         rotation = zeros(size(hamltMatrix)...)
         eigVals = zeros(size(hamltMatrix)[2])
         if !isempty(quantumNos)
-            for qantumNo in unique(quantumNos)
-                indices = findall(==(qantumNo), quantumNos)
+            numCaptured = 0
+            for quantumNo in unique(quantumNos)
+                indices = findall(==(quantumNo), quantumNos)
+                numCaptured += length(indices)
                 F = eigen(Hermitian(hamltMatrix[indices, indices]))
                 rotation[indices, indices] .= F.vectors
                 eigVals[indices] .= F.values
             end
+            @assert numCaptured == length(quantumNos)
             sortSequence = sortperm(eigVals)
             rotation = rotation[:, sortperm(eigVals)]
             quantumNos = quantumNos[sortperm(eigVals)]
@@ -146,6 +151,18 @@ function IterDiag(
                           for col in eachcol(rotation)]
         end
 
+        @time if step == length(hamltFlow)
+            serialize(savePaths[step], Dict("basis" => rotation,
+                                            "eigVals" => eigVals,
+                                            "quantumNos" => quantumNos,
+                                            "currentSites" => currentSites,
+                                            "newSites" => newSitesFlow[step],
+                                            "bondAntiSymmzer" => bondAntiSymmzer,
+                                           )
+                     )
+            break
+        end
+
         if isnothing(quantumNoReq) && length(eigVals) > maxSize
             # ensure we aren't truncating in the middle of degenerate states
             rotation = rotation[:, eigVals .≤ eigVals[maxSize] + abs(eigVals[maxSize]) * degenTol]
@@ -169,18 +186,6 @@ function IterDiag(
             eigVals = eigVals[retainBasisVectors]
         end
 
-        if step == length(hamltFlow)
-            serialize(savePaths[step], Dict("basis" => rotation,
-                                            "eigVals" => eigVals,
-                                            "quantumNos" => quantumNos,
-                                            "currentSites" => currentSites,
-                                            "newSites" => newSitesFlow[step],
-                                            "bondAntiSymmzer" => bondAntiSymmzer,
-                                           )
-                     )
-            break
-        end
-
         newBasis = BasisStates(length(newSitesFlow[step+1]))
 
         identityEnv = length(newSitesFlow[step+1]) == 1 ? I(2) : kron(fill(I(2), length(newSitesFlow[step+1]))...)
@@ -230,12 +235,12 @@ function IterDiag(
 
         # expand new operators
         bondAntiSymmzer = rotation' * bondAntiSymmzer * rotation
-        for site in newSitesFlow[step+1] 
+        @time for site in newSitesFlow[step+1] 
             operators[("+", [site])] = kron(bondAntiSymmzer, operators[("+", [site])])
             operators = CreateDNH(operators, site)
         end
 
-        operators = CreateProductOperator(create[step+1], operators, newSitesFlow[step+1])
+        @time operators = CreateProductOperator(create[step+1], operators, newSitesFlow[step+1])
 
         # expand the antisymmetrizer
         bondAntiSymmzer = kron(bondAntiSymmzer, fill(sigmaz, length(newSitesFlow[step+1]))...)