DBF Energy fluctuations

src/qibo/models/double_bracket.py

@@ -2,6 +2,8 @@
 
 import numpy as np
 
+from qibo.backends.numpy import NumpyBackend
+
 from ..config import raise_error
 from ..hamiltonians import Hamiltonian
 
@@ -82,16 +84,29 @@ def diagonal_h_matrix(self):
         """Diagonal H matrix."""
         return self.backend.cast(np.diag(np.diag(self.backend.to_numpy(self.h.matrix))))
 
+    @property
+    def off_diag_h(self):
+        return self.h.matrix - self.diagonal_h_matrix
+
     @property
     def off_diagonal_norm(self):
         """Norm of off-diagonal part of H matrix."""
-        off_diag_h = self.h.matrix - self.diagonal_h_matrix
         off_diag_h_dag = self.backend.cast(
-            np.matrix(self.backend.to_numpy(off_diag_h)).getH()
+            np.matrix(self.backend.to_numpy(self.off_diag_h)).getH()
+        )
+        return np.real(
+            np.trace(self.backend.to_numpy(off_diag_h_dag @ self.off_diag_h))
         )
-        return np.real(np.trace(self.backend.to_numpy(off_diag_h_dag @ off_diag_h)))
 
     @property
     def backend(self):
         """Get Hamiltonian's backend."""
         return self.h0.backend
+
+    def energy_fluctuation(self, state):
+        """Evaluate energy fluctuations"""
+        energy = self.h.expectation(state)
+        h = self.h.matrix
+        h2 = Hamiltonian(nqubits=self.h.nqubits, matrix=h @ h, backend=self.backend)
+        average_h2 = self.backend.calculate_expectation_state(h2, state, normalize=True)
+        return np.sqrt(average_h2 - energy**2)

tests/test_models_dbf.py

@@ -2,6 +2,7 @@
 import numpy as np
 import pytest
 
+from qibo.backends import GlobalBackend
 from qibo.hamiltonians import Hamiltonian
 from qibo.models.double_bracket import DoubleBracketFlow, FlowGeneratorType
 from qibo.quantum_info import random_hermitian
@@ -51,3 +52,11 @@ def test_double_bracket_flow_single_commutator(backend, nqubits):
         dbf(mode=FlowGeneratorType.single_commutator, step=0.01, d=d)
 
     assert initial_off_diagonal_norm > dbf.off_diagonal_norm
+
+
+def test_energy_fluctuations(backend):
+    h0 = np.array([[1, 0], [0, -1]])
+    state = np.array([1, 0])
+    dbf = DoubleBracketFlow(Hamiltonian(1, matrix=h0, backend=backend))
+    energy_fluctuation = dbf.energy_fluctuation(state=state)
+    assert energy_fluctuation == 0