MultiCouplingModel and Toric Code as example

Discussing the best way to implement feature X
Post Reply
User avatar
Johannes
Site Admin
Posts: 413
Joined: 21 Jul 2018, 12:52
Location: TU Munich

MultiCouplingModel and Toric Code as example

Post by Johannes »

Since I've recently added the MultiCouplingModel, it's time to add a simple model actually using it.
I think the exactly solvable toric code is a good example, so here is my implementation:

Code: Select all

"""Kitaev's exactly solvable toric code model.

As we put the model on a cylinder, the name "toric code" is a bit misleading,
but it is established for this model...
"""
# Copyright 2018 TeNPy Developers

import numpy as np

from .lattice import Lattice
from ..networks.site import SpinHalfSite
from .model import MultiCouplingModel, MPOModel
from ..tools.params import get_parameter, unused_parameters
from ..tools.misc import any_nonzero

__all__ = ['DualSquare', 'ToricCode']


class DualSquare(Lattice):
    """The dual lattice of the square lattice (again square).

    The sites in this lattice correspond to the vertical and horizontal (nearest neighbor) bonds
    of a common :class:`~tenpy.models.lattice.Square` lattice with the same dimensions `Lx, Ly`.

    Parameters
    ----------
    Lx, Ly : int
        Dimensions of the original lattice. This lattice has `2*Lx*Ly` sites.
    site_hor, site_ver : :class:`~tenpy.networks.site.Site`
        The sites for the horizontal and vertical bonds.
    order : str | (priority, snake_winding)
        A string or tuple specifying the order, given to :meth:`ordering`.
    """
    def __init__(self, Lx, Ly, site_hor, site_ver, order='default', bc_MPS='finite'):
        basis = np.eye(2)
        pos = np.array([[0.5, 0.], [0., 0.5]])
        self.nearest_neighbors = [(0, 1, np.array([0, 0])), (0, 1, np.array([1, 1])),
                                  (1, 0, np.array([-1, 1])), (1, 0, np.array([0, 1]))]
        self.next_nearest_neighbors = [(i, i, dx) for i in [0, 1]
                                       for dx in [np.array([1,0]), np.array([0, 1])]]
        super().__init__([Lx, Ly], [site_hor, site_ver], order, bc_MPS, basis, pos)


class ToricCode(MultiCouplingModel, MPOModel):
    r"""Spin-S sites coupled by nearest neighbour interactions.

    The Hamiltonian reads:

    .. math ::
        H = - \mathtt{Jv} \sum_{vertices v} \prod_{i \in v}  \sigma^x_i
            - \mathtt{Jp} \sum_{plaquettes p} \prod_{i \in p} \sigma^z_i

    (Note that this are Pauli matrices, not spin-1/2 operators.)
    All parameters are collected in a single dictionary `model_param` and read out with
    :func:`~tenpy.tools.params.get_parameter`.

    Parameters
    ----------
    Lx, Ly : int
        Dimension of the lattice, number of plaquettes around the cylinder.
    conserve : 'parity' | None
        What should be conserved. See :class:`~tenpy.networks.Site.SpinHalfSite`.
    Jc, Jp: float | array
        Couplings as defined for the Hamiltonian above.
    bc_MPS : {'finite' | 'infinte'}
        MPS boundary conditions. Coupling boundary conditions are chosen appropriately.
    order : str
        The order of the lattice sites in the lattice, see :class:`DualSquare`.
    """

    def __init__(self, model_param):
        # 0) read out/set default parameters
        verbose = get_parameter(model_param, 'verbose', 1, self.__class__)
        Lx = get_parameter(model_param, 'Lx', 2, self.__class__)
        Ly = get_parameter(model_param, 'Ly', 2, self.__class__)
        Jv = get_parameter(model_param, 'Jv', 1., self.__class__)
        Jp = get_parameter(model_param, 'Jp', 1., self.__class__)
        bc_MPS = get_parameter(model_param, 'bc_MPS', 'infinite', self.__class__)
        order = get_parameter(model_param, 'order', 'default', self.__class__)
        conserve = get_parameter(model_param, 'conserve', 'parity', self.__class__)
        unused_parameters(model_param, self.__class__)
        # 1) define Site and lattice
        site = SpinHalfSite(conserve)
        lat = DualSquare(Lx, Ly, site, site, order, bc_MPS)
        # 2) initialize CouplingModel
        bc_coupling = [None, 'periodic']
        bc_coupling[0] = 'periodic' if bc_MPS == 'infinite' else 'open'
        MultiCouplingModel.__init__(self, lat, bc_coupling)
        # 3) add terms of the Hamiltonian
        # (u is always 0 as we have only one site in the unit cell)
        Jv = np.asarray(Jv)
        Jp = np.asarray(Jp)
        # vertex/star term
        self.add_multi_coupling(Jv, 0, 'Sigmax', [(1, 'Sigmax', [0, 0]),
                                                  (0, 'Sigmax', [-1, 0]),
                                                  (1, 'Sigmax', [0, -1])])
        # plaquette term
        self.add_multi_coupling(Jp, 0, 'Sigmaz', [(1, 'Sigmaz', [0, 0]),
                                                  (0, 'Sigmaz', [0, 1]),
                                                  (1, 'Sigmaz', [1, 0])])
        # 4) initialize MPO
        MPOModel.__init__(self, lat, self.calc_H_MPO())

[mention]Umberto Borla[/mention] , can you do me the favor to check that the resulting energies are correct and that I didn't mess up any couplings?

Is this multi-coupling intuitive enough regarding how to use it? Any suggestions how to make it more intuitive?

Writing the code, I felt that it might be more intuitive to have a single list ops in add_multi_coupling, i.e. along the lines of

Code: Select all

self.add_multi_coupling(Jv, [(0, 'Sigmax', [0, 0]),
                             (1, 'Sigmax', [0, 0]),
                             (0, 'Sigmax', [-1, 0]),
                             (1, 'Sigmax', [0, -1])])
# instead of
self.add_multi_coupling(Jv, 0, 'Sigmax', [(1, 'Sigmax', [0, 0]),
                                          (0, 'Sigmax', [-1, 0]),
                                          (1, 'Sigmax', [0, -1])])
What do you think?
Umberto Borla
Posts: 18
Joined: 23 Jul 2018, 09:23
Location: Technical University Munich

Re: MultiCouplingModel and Toric Code as example

Post by Umberto Borla »

Hi [mention]Johannes[/mention] . I am looking at how the multi-coupling model works. However I have a quick question specific to the toric code example. What is L_x, when the boundary conditions are set to "infinite" (infinite cylinder)? I thought it would be just a repetition of the unit cell, and the size of the initial mps that I have to input is consistent with this assumption. However in that case I would imagine the energy per site not to depend on L_x, but it does. Am I missing something?
User avatar
Johannes
Site Admin
Posts: 413
Joined: 21 Jul 2018, 12:52
Location: TU Munich

Re: MultiCouplingModel and Toric Code as example

Post by Johannes »

You're right, L_x is just the number of "rings" of the cylinder included into the MPS unit cell. I agree that the energy per site shouldn't depend on L_x, so it's a bad sign if it does, there's probably something wrong then :(
User avatar
Johannes
Site Admin
Posts: 413
Joined: 21 Jul 2018, 12:52
Location: TU Munich

Re: MultiCouplingModel and Toric Code as example

Post by Johannes »

ok, which initial state and parameters are you using? Mixer enabled?
I checked with the following code

Code: Select all

import numpy as np
from tenpy.models.toric_code import ToricCode
from tenpy.networks.mps import MPS
from tenpy.algorithms import dmrg

for Lx in [1, 2, 3, 4]:
    print("="*80)
    print("Lx =", Lx)
    model_params = dict(Lx=Lx, verbose=0)
    M = ToricCode(model_params)
    psi = MPS.from_product_state(M.lat.mps_sites(), [0] * M.lat.N_sites, bc='infinite')

    dmrg_params = {
        'mixer': True,
        'trunc_params': {
            'chi_max': 50,
            'svd_min': 1.e-10
        },
        'max_E_err': 1.e-10,
        'verbose': 0
    }
    result = dmrg.run(psi, M, dmrg_params)
    print("E = {E:.13f}".format(E=result['E']))
    print("final bond dimensions: ", psi.chi)
and get (apart from some warnings) an energy E=-1 per site. Is that what you expect?
For the above parameters, I get different energies for some Lx and Ly=3, but I'm quite sure that's a DMRG convergence issue. The mixer seems to be critical, it might be better to enable the mixer for a few more sweeps with a bit more strength.

I guess you can write the exact state as an MPS, for Ly=2 the bond dimensions are just [2, 4, 4, 4].
The eigenvalues of the transfer matrix are 1 and 0 up to rounding errors, is that correct?
Could you check the overlap?
Umberto Borla
Posts: 18
Joined: 23 Jul 2018, 09:23
Location: Technical University Munich

Re: MultiCouplingModel and Toric Code as example

Post by Umberto Borla »

I tested your code a bit more and compared to mine, setting \(L_x = 1\) in both cases (as we agree, it should not matter). What I found out is that the results coincide exactly with mine whenever \(L_y\) is even, and fails otherwise (e.g. the energy depends on \(L_x\), the bond dimension is always saturated, entanglement entropy is incorrect and all sort of bad things...). So I guess it must be just a strange bug. Any ideas? Also, I observed that your code is much faster due to parity conservation.

I think that my code is correct because it reproduces the following features of the TC for any \(L_y > 2\):

- Ground state energy (it is easily calculated analytically, we just have to agree if we want the energy per site or per cell)
- Entanglement entropy: it is \((L-1)\ln 2\), the minus one being the topological contribution.
- The maximum bond dimension of the MPS is \(2^{L_y + 1}\)
- The correlation length is very close to zero

Just to be sure, this is how I run your code:

Code: Select all

def example_run_TC():
    model_params = dict(Lx=1, Ly=3, Jv=1, Jp=1, bc_MPS="infinite", verbose= 1)
    Lx = model_params['Lx']
    Ly = model_params['Ly']
    model = ToricCode(model_params)
    psi = MPS.from_product_state(model.lat.mps_sites(), [0]*(2*Ly*Lx), bc = "infinite")
    dmrg_params = {'mixer': True,
                   'chi_list': {0: 80},
                   'trunc_params': {'svd_min': 1.e-10},
                   'verbose': 1}
    results = dmrg.run(psi, model, dmrg_params)
    print("Energy per site: ", results['E'])
    print("Entanglement Entropy: ", psi.entanglement_entropy()[0])

if __name__ == "__main__":
    example_run_TC()
Umberto Borla
Posts: 18
Joined: 23 Jul 2018, 09:23
Location: Technical University Munich

Re: MultiCouplingModel and Toric Code as example

Post by Umberto Borla »

Oh I just saw your previous reply. Indeed, for L_y = 2 we get the same answers. The problem seems to be when L_y is odd.
User avatar
Johannes
Site Admin
Posts: 413
Joined: 21 Jul 2018, 12:52
Location: TU Munich

Re: MultiCouplingModel and Toric Code as example

Post by Johannes »

I checked explicitly, that the Hamiltonian is correct for Lx=1, Ly=3 by multipliying out the MPO matrices. I get the terms
\[X_0 X_6 X_7 X_{11} + X_2 X_7 X_8 X_9 + X_4 X_9 X_{10} X_{11} + Z_0 Z_1 Z_2 Z_7 + Z_2 Z_3 Z_4 Z_9 + Z_0 Z_4 Z_5 Z_{11}\]
starting in the MPS unit cell of sites 0-5, which is exactly what I expected.
In other words, it really is a convergence problem of DMRG being stuck in a local minimum.
It might be solved by choosing another initial state. Since there is an analytical solution, it's porbably also a good idea to write a function initializing the exact ground state as an MPS (did you write something like that already for youself?). It can also serve as a good guess when you perturb the toric code model with additional terms.
In the long-term, another solution might be to enhance the "Mixer" of DMRG, this might also be useful for other (2D) models. It's on the TODO-list...
Umberto Borla
Posts: 18
Joined: 23 Jul 2018, 09:23
Location: Technical University Munich

Re: MultiCouplingModel and Toric Code as example

Post by Umberto Borla »

I did some tests myself and i reached the same conclusion, that the MPO matrices give the correct Hamiltonian once multiplied together (I also tried explicitly for Ly = 3).

I used the MPO matrices of TeNpy's toric-code model to create another toric-code model using the grid_outer method. For Lx=1, Ly=3, this gives the same convergence problems as the TeNpy model, as expected. However for Lx=2, Ly=3 this gives the correct value 2*log2 for the entanglement entropy, while the TeNpy model despite converging to the correct energy, gives 3*log2 for the entropy. This is also reflected by the maximum MPS bond dimension, which is 16 in the first case and 32 in the second. Does this give you any ideas?

Note that my original MPO matrices, which are unoptimized and give "disconnected" MPOs as you noted, lead to the correct result for the entropy and DMRG converges also for Lx=1. Could it be that there is some subtle loss of information in TeNpy's routine that creates optimized MPOs?
Post Reply