From 8cda315f1f9a78fc3045101ca827cd4b4c2a6a70 Mon Sep 17 00:00:00 2001
From: Quantumgunhee <gunheepark0229@gmail.com>
Date: Thu, 11 Jan 2024 11:59:51 -0800
Subject: [PATCH] mps get bipartite ent

---
 pyblock3/algebra/mps.py | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/pyblock3/algebra/mps.py b/pyblock3/algebra/mps.py
index 02dfe1e..38f2145 100644
--- a/pyblock3/algebra/mps.py
+++ b/pyblock3/algebra/mps.py
@@ -366,6 +366,25 @@ def compress(self, **opts):
                 tensors[i + 1] = np.tensordot(sr, tensors[i + 1], axes=1)
                 merror = max(merror, err)
         return type(self)(tensors=tensors, const=self.const, opts=self.opts, dq=self.dq), merror
+    
+    def get_bipartite_entanglement(self):
+        tensors = [ts for ts in self.tensors]
+        bip_ent = np.zeros(self.n_sites - 1, dtype=np.float64)
+        for i in range(self.n_sites - 1, 0, -1):
+            l, q = tensors[i].right_canonicalize()
+            tensors[i] = q
+            tensors[i - 1] = np.tensordot(tensors[i - 1], l, axes=1)
+        for i in range(0, self.n_sites - 1, 1):
+            l, s, r = tensors[i].left_svd(full_matrices=False)
+            sr = np.tensordot(s.diag(), r, axes=1)
+            tensors[i] = l
+            tensors[i + 1] = np.tensordot(sr, tensors[i + 1], axes=1)
+
+            ldsq = s.data**2
+            ldsq /= np.sum(ldsq)
+            ldsq = ldsq[ldsq != 0]
+            bip_ent[i] = float(np.sum(-ldsq * np.log(ldsq)))
+        return bip_ent
 
     @classmethod
     def _add(cls, a, b):