Merge pull request #171 from UCL-CCS/state_rdm

State rdm

symmer/operators/base.py

@@ -1068,8 +1068,9 @@ def commutes(self,
         Returns:
             bool: True if all terms commute, False otherwise.
         """
-        return self.commutator(PwordOp).n_terms == 0
-
+        commutator = self.commutator(PwordOp).cleanup()
+        return (commutator.n_terms == 0 or np.all(commutator.coeff_vec[0] == 0))
+
     @cached_property
     def adjacency_matrix(self) -> np.array:
         """ 
@@ -2016,6 +2017,45 @@ def to_sparse_matrix(self):
             # conjugate has already taken place, just need to make into row vector
             sparse_Qstate= sparse_Qstate.reshape([1,-1])
         return sparse_Qstate
+
+    @cached_property
+    def to_dense_matrix(self):
+        """
+        Returns:
+            dense_Qstate (ndarray): dense matrix representation of the statevector
+        """
+        return self.to_sparse_matrix.toarray()
+
+    def partial_trace_over_qubits(self, qubits: List[int] = []) -> np.ndarray:
+        """
+        Perform a partial trace over the specified qubit positions, 
+        yielding the reduced density matrix of the remaining subsystem.
+
+        Args:
+            qubits (List[int]): qubit indicies to trace over
+
+        Returns:
+            rho_reduced (ndarray): Reduced density matrix over the remaining subsystem
+        """
+        rho_reduced = self.to_dense_matrix.reshape([2]*self.n_qubits)
+        rho_reduced = np.tensordot(rho_reduced, rho_reduced.conj(), axes=(qubits, qubits))
+        d = int(np.sqrt(np.product(rho_reduced.shape)))
+        return rho_reduced.reshape(d, d)
+
+    def get_rdm(self, qubits: List[int] = []) -> np.ndarray:
+        """
+        Return the reduced density matrix of the specified qubit positions, 
+        corresponding with a partial trace over the complementary qubit indices
+        
+        Args:
+            qubits (List[int]): qubit indicies to preserve
+
+        Returns:
+            rho_reduced (ndarray): Reduced density matrix over the chosen subsystem
+        """
+        trace_over_indices = list(set(range(self.n_qubits)).difference(set(qubits)))
+        rho_reduced = self.partial_trace_over_qubits(trace_over_indices)
+        return rho_reduced
 
     def _is_normalized(self) -> bool:
         """

symmer/operators/noncontextual_op.py

@@ -551,12 +551,12 @@ def get_symmetry_contributions(self, nu: np.array) -> float:
         si = np.array([np.sum(coeff_mod[mask]).real for mask in self.mask_Ci])
         return s0, si
 
-    def get_energy(self, nu: np.array) -> float:
-        """
+    def get_energy(self, nu: np.array, AC_ev: int = -1) -> float:
+        """ 
         The classical objective function that encodes the noncontextual energies.
         """
         s0, si = self.get_symmetry_contributions(nu)
-        return s0 - np.linalg.norm(si)
+        return s0 + AC_ev * np.linalg.norm(si, ord=2)
 
     def update_clique_representative_operator(self, clique_index:int = None) -> List[Tuple[PauliwordOp, float]]:
         _, si = self.get_symmetry_contributions(self.symmetry_generators.coeff_vec)

symmer/utils.py

@@ -75,6 +75,23 @@ def exact_gs_energy(
         # if a solution is not found within the first n_eig eigenvalues then error
         raise RuntimeError('No eigenvector of the correct particle number was identified - try increasing n_eigs.')
 
+def get_entanglement_entropy(psi: QuantumState, qubits: List[int]) -> float:
+    """
+    Get the Von Neumann entropy of the biprtition defined by the specified subsystem 
+    qubit indices and those remaining (i.e. those that will be subsequently traced out)
+
+    Args:
+        psi (QuantumState): the quantum state for which we wish to extract the entanglement entropy
+        qubits (List[int]): the qubit indices to project onto (the remaining qubits will be traced over)
+    
+    Returns:
+        entropy (float): the Von Neumann entropy of the reduced subsystem
+    """
+    reduced = psi.get_rdm(qubits)
+    eigvals, eigvecs = np.linalg.eig(reduced)
+    eigvals = eigvals[eigvals>0]
+    entropy = -np.sum(eigvals*np.log(eigvals)).real
+    return entropy
 
 def random_anitcomm_2n_1_PauliwordOp(n_qubits, complex_coeff=False, apply_clifford=True):
     """