implement PurificationMPS.from_infiniteT_canonical()

doc/changelog/latest.rst

@@ -45,6 +45,7 @@ Added
 - :attr:`tenpy.models.lattice.Lattice.Lu` as a class attribute.
 - :meth:`tenpy.models.lattice.Lattice.find_coupling_pairs` to automatically find coupling pairs of 'nearest_neighbors' etc..
 - :class:`tenpy.models.lattice.HelicalLattice` allowing to have a much smaller MPS unit cell by shifting the boundary conditions around the cylinder.
+- :meth:`tenpy.networks.purification_mps.PurificationMPS.from_infiniteT_canonical` for a canonical ensemble.
 
 Changed
 ^^^^^^^

tenpy/algorithms/exact_diag.py

@@ -237,7 +237,7 @@ def mps_to_full(self, mps):
         Returns
         -------
         psi : :class:`~tenpy.linalg.np_conserved.Array`
-            The MPO contracted along the virtual bonds.
+            The MPS contracted along the virtual bonds.
         """
         if mps.bc != 'finite':
             raise ValueError("Full diagonalization works only on finite systems")

tenpy/networks/purification_mps.py

@@ -102,17 +102,11 @@
     conservation with the above `qconj` convention.
     Moreover, we don't split the physical and auxiliar space into separate sites, which makes
     TEBD as costly as :math:`O(d^6 \chi^3)`.
-
-.. todo ::
-    One can also look at the canonical ensembles by defining the conserved quantities
-    differently, see [barthel2016]_ for details.
-    Idea: usual charges on `p`, trivial charges on `q`; fix total charge to desired value.
-    I think it should suffice to implement another `from_infiniteT`.
-
 """
 # Copyright 2018-2021 TeNPy Developers, GNU GPLv3
 
 import numpy as np
+import itertools
 
 from .mps import MPS
 from ..linalg import np_conserved as npc
@@ -154,18 +148,21 @@ def test_sanity(self):
         super().test_sanity()
 
     @classmethod
-    def from_infiniteT(cls, sites, bc='finite', form='B'):
-        """Initial state corresponding to infinite-Temperature ensemble.
+    def from_infiniteT(cls, sites, bc='finite', form='B', dtype=np.float64):
+        """Initial state corresponding to grand-canonical infinite-temperature ensemble.
 
         Parameters
         ----------
         sites : list of :class:`~tenpy.networks.site.Site`
             The sites defining the local Hilbert space.
+            For usual :class:`tenpy.models.model.Model` given by `model.lat.mps_sites()`.
         bc : {'finite', 'segment', 'infinite'}
             MPS boundary conditions as described in :class:`~tenpy.networks.mps.MPS`.
         form : (list of) {``'B' | 'A' | 'C' | 'G' | None`` | tuple(float, float)}
             The canonical form of the stored 'matrices', see table in :mod:`~tenpy.networks.mps`.
             A single choice holds for all of the entries.
+        dtype : type or string
+            The data type of the array entries.
 
         Returns
         -------
@@ -179,13 +176,100 @@ def from_infiniteT(cls, sites, bc='finite', form='B'):
         Bs = [None] * L
         for i in range(L):
             p_leg = sites[i].leg
-            B = npc.diag(1., p_leg, np.float, ['p', 'q']) / sites[i].dim**0.5
+            B = npc.diag(1., p_leg, dtype, ['p', 'q']) / sites[i].dim**0.5
             # leg `q` has the physical leg with opposite `qconj`
             B = B.add_trivial_leg(0, label='vL', qconj=+1).add_trivial_leg(1, label='vR', qconj=-1)
             Bs[i] = B
         res = cls(sites, Bs, S, bc, form)
         return res
 
+    @classmethod
+    def from_infiniteT_canonical(cls, sites, charge_sector, form='B', dtype=np.float64):
+        """Initial state corresponding to *canonical* infinite-temperature ensemble.
+
+        Works only for finite boundary conditions, following the idea outlined in [barthel2016]_.
+        However, we just put trivial charges on the ancilla legs,
+        instead of doubling the number of charges as suggested in that paper.
+
+        Note that at least some of the disentanglers don't work with the canonical ensemble.
+
+        Parameters
+        ----------
+        sites : list of :class:`~tenpy.networks.site.Site`
+            The sites defining the local Hilbert space.
+            For usual :class:`tenpy.models.model.Model` given by `model.lat.mps_sites()`.
+        charge_sector : tuple of int
+            The desired charge sector to be taken for the canonical ensemble.
+        form : (list of) {``'B' | 'A' | 'C' | 'G' | None`` | tuple(float, float)}
+            The canonical form of the stored 'matrices', see table in :mod:`~tenpy.networks.mps`.
+            A single choice holds for all of the entries.
+
+        Returns
+        -------
+        infiniteT_MPS : :class:`PurificationMPS`
+            Describes the infinite-temperature (grand canonical) ensemble,
+            i.e. expectation values give a trace over all basis states.
+        """
+        sites = list(sites)
+        L = len(sites)
+        assert L > 0
+        chinfo = sites[0].leg.chinfo
+        for s in sites:
+            assert s.leg.chinfo == chinfo
+        charge_sector_left = chinfo.make_valid(None)  # zero charges
+        charge_sector_right = chinfo.make_valid(charge_sector)
+        assert charge_sector_right.ndim == 1
+        # get bounds for the maximal and minimal charge values at each bond
+        qflats = [s.leg.to_qflat() for s in sites]
+        min_p_Q = [np.min(qflat, axis=0) for qflat in qflats]
+        max_p_Q = [np.max(qflat, axis=0) for qflat in qflats]
+        min_Q_left = charge_sector_left + np.cumsum(min_p_Q, axis=0)  # on bonds 1, 2, 3... L
+        max_Q_left = charge_sector_left + np.cumsum(max_p_Q, axis=0)
+        min_Q_right = charge_sector_right - np.cumsum(max_p_Q[::-1], axis=0)[::-1]  # 0, 1, ... L-1
+        max_Q_right = charge_sector_right - np.cumsum(min_p_Q[::-1], axis=0)[::-1]
+        min_Q = np.max([min_Q_left[:-1], min_Q_right[1:]], axis=0)  # on non-trivial bonds
+        max_Q = np.min([max_Q_left[:-1], max_Q_right[1:]], axis=0)  # on non-trivial bonds
+        min_Q = np.append(min_Q, [charge_sector_right], axis=0)  # bonds 1, 2, ... L
+        max_Q = np.append(max_Q, [charge_sector_right], axis=0)
+        assert np.all(max_Q >= min_Q)
+        chi = np.prod(max_Q - min_Q + 1, axis=1)  # assumes that charges can be incremented by 1
+
+        Bs = []
+        Ss = [np.ones(1, dtype=np.float64)]
+        # now we can define the tensors following [barthel2016]_
+        # B[vL, vR, p, q] = delta_{p,q} delta_{Q(p) + Q(vL), Q(vR)}
+        right_Q = chinfo.make_valid([charge_sector_left])
+        right_leg = npc.LegCharge.from_qflat(chinfo, right_Q, qconj=-1)
+        for s in range(L):
+            p_leg = sites[s].leg
+            p_Q = p_leg.to_qflat()
+            q_leg = p_leg.conj().copy()
+            q_leg.charges = np.zeros_like(q_leg.charges)
+            left_Q = right_Q
+            left_leg = right_leg.conj()
+            right_Q = np.array(
+                list(
+                    itertools.product(
+                        *[range(min_Q[s][c], max_Q[s][c] + 1) for c in range(chinfo.qnumber)])))
+            right_Q = right_Q[np.lexsort(right_Q.T), :]  # sort charges
+            right_leg = npc.LegCharge.from_qflat(chinfo, right_Q, qconj=-1)
+            B = npc.zeros([left_leg, right_leg, p_leg, q_leg],
+                          dtype=dtype,
+                          labels=['vL', 'vR', 'p', 'q'])
+            for p in range(p_leg.ind_len):
+                for vL in range(left_Q.shape[0]):
+                    Q_vR = left_Q[vL] + p_Q[p]
+                    vR = np.nonzero(np.all(Q_vR == right_Q, axis=1))[0]
+                    if len(vR) == 0:
+                        continue
+                    vR = vR.item()
+                    B[vL, vR, p, p] = 1.  # add an entry in the tensor
+            Bs.append(B)
+            Ss.append(np.ones(B.shape[1], np.float64))
+        res = cls(sites, Bs, Ss, 'finite', form)
+        res.canonical_form_finite()  # calculate S values
+        return res
+
     def entanglement_entropy_segment(self, segment=[0], first_site=None, n=1, legs='p'):
         r"""Calculate entanglement entropy for general geometry of the bipartition.
 

tests/test_purification.py

@@ -4,6 +4,7 @@
 import warnings
 import numpy as np
 import numpy.testing as npt
+import scipy
 from tenpy.models.xxz_chain import XXZChain
 
 from tenpy.networks import purification_mps, site
@@ -14,6 +15,7 @@
 import pytest
 
 spin_half = site.SpinHalfSite(conserve='Sz')
+ferm = site.FermionSite(conserve='N')
 
 
 def test_purification_mps():
@@ -44,6 +46,47 @@ def test_purification_mps():
             npt.assert_allclose(C, 0.5 * 0.5 * np.eye(L), atol=1.e-13)
 
 
+def test_canoncial_purification(L=6, charge_sector=0, eps=1.e-14):
+    site = spin_half
+    psi = purification_mps.PurificationMPS.from_infiniteT_canonical([site] * L, [charge_sector])
+    psi.test_sanity()
+    total_psi = psi.get_theta(0, L).take_slice(0, 'vL').take_slice(0, 'vR')
+    total_psi.itranspose(['p' + str(i) for i in range(L)] + ['q' + str(i) for i in range(L)])
+    # note: don't `combine_legs`: it will permute the p legs differently than q due to charges
+    total_psi_dense = total_psi.to_ndarray().reshape(2**L, 2**L)
+    # now it should be diagonal
+    diag = np.diag(total_psi_dense)
+    assert np.all(np.abs(total_psi_dense - np.diag(diag) < eps))  # is it diagonal?
+    # and the diagonal should be sqrt(L choose L//2) for states with fitting numbers
+    pref = 1. / scipy.special.comb(L, L // 2 + charge_sector)**0.5
+    Q_p = site.leg.to_qflat()[:, 0]
+    for i, entry in enumerate(diag):
+        Q_i = sum([Q_p[int(b)] for b in format(i, 'b').zfill(L)])
+        if Q_i == charge_sector:
+            assert abs(entry - pref) < eps
+        else:
+            assert abs(entry) < eps
+
+    # and one quick test of TEBD
+    xxz_pars = dict(L=L, Jxx=1., Jz=3., hz=0., bc_MPS='finite')
+    M = XXZChain(xxz_pars)
+    TEBD_params = {
+        'trunc_params': {
+            'chi_max': 16,
+            'svd_min': 1.e-8
+        },
+        'disentangle': None,  # 'renyi' should work as well, 'backwards' not.
+        'dt': 0.1,
+        'verbose': 30,
+        'N_steps': 2
+    }
+    eng = PurificationTEBD(psi, M, TEBD_params)
+    eng.run_imaginary(0.2)
+    eng.run()
+    N = psi.expectation_value('Id')  # check normalization : <1> =?= 1
+    npt.assert_array_almost_equal_nulp(N, np.ones([L]), 100)
+
+
 @pytest.mark.slow
 def test_purification_TEBD(L=3):
     xxz_pars = dict(L=L, Jxx=1., Jz=3., hz=0., bc_MPS='finite')