implement MPS.apply_product_op()

doc/changelog/latest.rst

@@ -46,6 +46,7 @@ Added
 - Simulation class :class:`~tenpy.simulation.simulation.Simulation` and subclasses as a new extra layer for handling the general setup.
 - Command line script ``tenpy-run`` and :func:`~tenpy.run_simulation` for setting up a simulation.
 - :meth:`~tenpy.networks.mps.MPS.entanglement_entropy_segment2`
+- :meth:`~tenpy.networks.mps.MPS.apply_product_op`
 - :meth:`tenpy.linalg.sparse.FlatLinearOperator.eigenvectors` and :meth:`~tenpy.linalg.sparse.FlatHermitianOperator.eigenvectors` to unify
   code from :meth:`tenpy.networks.mps.TransferMatrix.eigenvectors` and :meth:`tenpy.linalg.lanczos.lanczos_arpack`.
 - :meth:`tenpy.tools.misc.group_by_degeneracy`

tenpy/networks/mps.py

@@ -2920,6 +2920,51 @@ def apply_local_op(self, i, op, unitary=None, renormalize=False, cutoff=1.e-13):
         if not unitary:
             self.canonical_form(renormalize)
 
+    def apply_product_op(self, ops, unitary=None, renormalize=False):
+        """Apply a (global) product of local onsite operators to `self`.
+
+        Note that this destroys the canonical form if any local operator is non-unitary.
+        Therefore, this function calls :meth:`canonical_form` if necessary.
+
+        The result is equivalent to the following loop, but more efficient by avoiding
+        intermediate calls to :meth:`canonical_form` inside the loop::
+
+            for i, op in enumerate(ops):
+                self.apply_local_op(i, op, unitary, renormalize, cutoff)
+
+        Parameters
+        ----------
+        ops : (list of) str | npc.Array
+            List of onsite operators to apply on each site, with legs ``'p', 'p*'``.
+            Strings (like ``'Id', 'Sz'``) are translated into single-site operators defined by
+            :attr:`sites`.
+        unitary : None | bool
+            Whether `op` is unitary, i.e., whether the canonical form is preserved (``True``)
+            or whether we should call :meth:`canonical_form` (``False``).
+            ``None`` checks whether ``max(norm(op dagger(op) - identity) for op in ops) < 1.e-14``
+        renormalize : bool
+            Whether the final state should keep track of the norm (False, default) or be
+            renormalized to have norm 1 (True).
+        """
+        ops = to_iterable(ops)
+        if self.L % len(ops) != 0:
+            raise ValueError("len of ops incommensurate with self.L")
+        self.convert_form('B')
+        for i in range(self.L):
+            op = ops[i % len(ops)]
+            if isinstance(op, str):
+                if op == 'Id':
+                    continue  # nothing to do here...
+                op = self.sites[i].get_op(op)
+            if unitary is None:
+                op_op_dagger = npc.tensordot(op, op.conj(), axes=['p*', 'p'])
+                if npc.norm(op_op_dagger - npc.eye_like(op_op_dagger)) > 1.e-14:
+                    unitary = False
+            # actually apply the operator at site i
+            self._B[i] = npc.tensordot(op, self._B[i], axes=['p*', 'p'])
+        if not unitary:
+            self.canonical_form(renormalize)
+
     def swap_sites(self, i, swap_op='auto', trunc_par=None):
         r"""Swap the two neighboring sites `i` and `i+1` (inplace).
 

tests/test_mps.py

@@ -271,6 +271,26 @@ def test_canonical_form(bc):
     assert np.max(psi.norm_test()) < 1.e-14
 
 
+@pytest.mark.parametrize("bc", ['finite', 'infinite'])
+def test_apply_op(bc, eps=1.e-13):
+    s = site.SpinHalfSite(None)
+    psi0 = mps.MPS.from_singlets(s, 3, [(0, 2)], lonely=[1], bc=bc, lonely_state='up')
+    psi1 = psi0.copy()
+    psi1.apply_local_op(1, 'Sigmax')  #unitary
+    psi1_expect = mps.MPS.from_singlets(s, 3, [(0, 2)], lonely=[1], bc=bc, lonely_state='down')
+    psi1 = psi0.copy()
+    psi1.apply_local_op(1, 'Sm')  #non-unitary
+    assert abs(psi1_expect.overlap(psi1) - 1.) < eps
+
+    psi2 = psi0.copy()
+    th = psi2.get_theta(0, 3).to_ndarray().reshape((8, ))
+    s2 = 0.5**0.5
+    assert np.linalg.norm(th - [0., s2, 0., 0., -s2, 0., 0, 0.]) < eps
+    psi2.apply_product_op(['Sigmax', 'Sm', 'Sigmax'])
+    th = psi2.get_theta(0, 3).to_ndarray().reshape((8, ))
+    assert np.linalg.norm(th - [0., 0., 0., -s2, 0., 0., s2, 0.]) < eps
+
+
 def test_enlarge_mps_unit_cell():
     s = site.SpinHalfSite(conserve='Sz')
     psi = mps.MPS.from_product_state([s] * 3, ['up', 'down', 'up'], bc='infinite')