fix complex Hamiltonians

Author: hczhai

docs/source/tutorial/mpo-mps-quimb.ipynb

@@ -213,7 +213,7 @@
     "    counter[0] += 1\n",
     "    return pympo.tensors[counter[0] - 1].blocks[0].reduced\n",
     "\n",
-    "qumpo = qtn.MatrixProductOperator.from_fill_fn(mpo_fill_fn, len(pyket.tensors), bond_dim=None, phys_dim=2, shape=\"ludr\")\n",
+    "qumpo = qtn.MatrixProductOperator.from_fill_fn(mpo_fill_fn, len(pyket.tensors), bond_dim=None, phys_dim=2, shape=\"ldur\")\n",
     "qumps = qtn.MatrixProductState.from_fill_fn(mps_fill_fn, len(pympo.tensors), bond_dim=None, phys_dim=2, shape=\"lpr\")\n",
     "\n",
     "qumps.show()\n",

pyblock2/algebra/io.py

@@ -251,6 +251,11 @@ def from_block2(bmps):
                 bmps.unload_tensor(i)
                 mr.deallocate()
                 r.deallocate()
+                if i == bmps.n_sites - 1 and tensors[i].rank == 3:
+                    for block in tensors[i].blocks:
+                        block.reduced = block.reduced.reshape(block.reduced.shape[:2])
+                        block.q_labels = block.q_labels[:2]
+                        block.rank = 2
             elif (
                 i == bmps.center and i != 0 and i != bmps.n_sites - 2 and bmps.dot == 2
             ):

pyblock2/driver/core.py

@@ -4471,6 +4471,7 @@ def dmrg(
         fused_contraction_multiplication=False,
         fused_contraction_rotation=False,
         lowmem_numerical_transform=False,
+        real_density_matrix=False,
     ):
         """
         Perform the ground state and/or excited state Density Matrix
@@ -4604,6 +4605,9 @@ def dmrg(
             lowmem_numerical_transform : bool
                 If True, numerical transform will be done with copying to save peak memory.
                 Defult is False.
+            real_density_matrix : bool
+                If True, will use the real part of the density matrix for truncation.
+                Should not be used for complex Hamiltonians. Default is False.
 
         Returns:
             energy : float|complex or list[float|complex]
@@ -4729,7 +4733,8 @@ def dmrg(
         dmrg.store_seq_data = store_seq_data
         dmrg.iprint = iprint
         dmrg.cutoff = cutoff
-        dmrg.trunc_type = dmrg.trunc_type | bw.b.TruncationTypes.RealDensityMatrix
+        if real_density_matrix:
+            dmrg.trunc_type = dmrg.trunc_type | bw.b.TruncationTypes.RealDensityMatrix
         if kernel is not None:
             # need to keep Python derived class in memory (stored in self)
             self.dmrg_kernel = self.make_kernel(kernel=kernel)