diff --git a/double_bracket_flow_as_a_diagonalization_quantum_algorithm b/double_bracket_flow_as_a_diagonalization_quantum_algorithm
new file mode 160000
index 0000000000..808dd6325b
--- /dev/null
+++ b/double_bracket_flow_as_a_diagonalization_quantum_algorithm
@@ -0,0 +1 @@
+Subproject commit 808dd6325b327a756a4b6ae5a4d560e01d1790cb
diff --git a/src/qibo/backends/numpy.py b/src/qibo/backends/numpy.py
index 62d4f5f5f5..7b9488b8fc 100644
--- a/src/qibo/backends/numpy.py
+++ b/src/qibo/backends/numpy.py
@@ -758,15 +758,15 @@ def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None):
         return self.np.matmul(eigenvectors, self.np.matmul(expd, ud))
 
     def calculate_expectation_state(self, hamiltonian, state, normalize):
-        ev = self.np.real(self.np.sum(self.np.conj(state) * hamiltonian.matrix @ state))
+        statec = self.np.conj(state)
+        hstate = hamiltonian @ state
+        ev = self.np.real(self.np.sum(statec * hstate))
         if normalize:
             ev /= self.np.sum(self.np.square(self.np.abs(state)))
         return ev
 
     def calculate_expectation_density_matrix(self, hamiltonian, state, normalize):
-        statec = self.np.conj(state)
-        hstate = hamiltonian @ state
-        ev = self.np.real(self.np.sum(statec * hstate))
+        ev = self.np.real(self.np.trace(self.cast(hamiltonian @ state)))
         if normalize:
             norm = self.np.real(self.np.trace(state))
             ev /= norm