Lattice design issue

[mrader]

This is not really a feature request, but I didn't know where else I should put this.

I wanted to bring up that the Lattice class has a severe design issue: It contains information about sites (which I would consider as information about the model) and also information about symmetries that will be used (which I would consider as information about the MPS). Would it be possible to separate these things? If you could do this, I guess it would not be neccessary anymore to implement the same model again and again for each lattice, but it would be possible to do for instance the following:

lattice = SquareLattice(Lx=3, Ly=4, bc_x="infinite", bc_y="cylinder")
model = TFIModel(lattice, J=-1, g=0.5)

[Johannes]

On one hand I agree that this would be very cool to have an single "TFIModel" with some specified nearest neighbor couplings.

On the other hand, I see some problems with that:

1. When you only want to derive the Model from the MPOModel, this approach is fine, but what if you want to derive from the NearestNeighborModel?
   Whether the model is actually nearest-neighbor (in the MPS sense) depends on the chosen lattice.
2. Where would you specify the sites in this approach? I would say that a 2-site ladder of two coupled spin-1/2 sites is a different geometry than a ladder coupling a spin-1/2 site with a spin-1 site; that's why I feel that the sites belong to the lattice class. I guess that's a convention what's geometry, though.
3. Once you have more than trivial nearest-neighbor couplings, I think it's better to have a fixed geometry in mind, especially when you need to think about the "dx" and "u" parameters for the "add_coupling" or when you want to create site-dependent couplings.

Given that the CouplingModel makes it so easy to create new models, I vote to rather keep it explicit with one lattice per model for most cases.
For simple cases, one can add a model parameter "lattice", though, as in the following example:

Python: Select all

class TFIModel2D(CouplingModel, MPOModel):
    r"""Transverse field Ising model on a 2D lattice.

    The Hamiltonian reads:

    .. math ::
        H = - \sum_{\langle i,j\rangle, i < j} \mathtt{J} \sigma^x_i \sigma^x_{j}
            - \sum_{i} \mathtt{g} \sigma^z_i

    Here, :math:`\langle i,j \rangle, i< j` denotes nearest neighbor pairs, each pair appearing
    exactly once.
    All parameters are collected in a single dictionary `model_param` and read out with
    :func:`~tenpy.tools.params.get_parameter`.

    Parameters
    ----------
    lattice : str | :class:`~tenpy.models.lattice.Lattice`
        The lattice class for the underlaying geometry. A string should be one of the lattices
        defined in th :mod:`~tenpy.models.lattice`.
    Lx, Ly : int
        Length of the lattice in x- and y-direction.
    J, g : float | array
        Couplings as defined for the Hamiltonian above.
    bc_MPS : {'finite' | 'infinte'}
        MPS boundary conditions along the x-direction.
        For 'infinite' boundary conditions, repeat the unit cell in x-direction.
        Coupling boundary conditions in x-direction are chosen accordingly.
    bc_y : 'ladder' | 'cylinder'
        Boundary conditions in y-direction.
    conserve : None | 'parity'
        What should be conserved. See :class:`~tenpy.networks.Site.SpinSite`.
    order : string
        Ordering of the sites in the MPS, e.g. 'default', 'snake';
        see :meth:`~tenpy.models.lattice.Lattice.ordering`.
    """

    def __init__(self, model_param):
        # 0) read out/set default parameters
        lat = get_parameter(model_param, 'lattice', "Square", self.__class__)
        Lx = get_parameter(model_param, 'Lx', 1, self.__class__)
        Ly = get_parameter(model_param, 'Ly', 4, self.__class__)
        J = get_parameter(model_param, 'J', 1., self.__class__)
        g = get_parameter(model_param, 'g', 1., self.__class__)
        bc_MPS = get_parameter(model_param, 'bc_MPS', 'infinite', self.__class__)
        bc_y = get_parameter(model_param, 'bc_y', 'cylinder', self.__class__)
        order = get_parameter(model_param, 'order', 'default', self.__class__)
        conserve = get_parameter(model_param, 'conserve', None, self.__class__)
        assert conserve != 'Sz'  # invalid!
        assert bc_y in ['cylinder', 'ladder']
        unused_parameters(model_param, self.__class__)  # checks for mistyped parameters
        # 1-3)
        site = SpinHalfSite(conserve=conserve)
        # 4) lattice
        if isinstance(lat, str):
            lat = lattice.__dict__.get(lat)  # the corresponding class
            lat = lat(Lx, Ly, site, order=order, bc_MPS=bc_MPS)  # instance of the class
        else:
            assert lat.Ls == (Lx, Ly)
        bc_coupling_x = 'periodic' if bc_MPS == 'infinite' else 'open'
        bc_coupling_y = 'periodic' if bc_y == 'cylinder' else 'open'
        # 5) initialize CouplingModel
        CouplingModel.__init__(self, lat, [bc_coupling_x, bc_coupling_y])
        # 6) add terms of the Hamiltonian
        # (u is always 0 as we have only one site in the unit cell)
        self.add_onsite(-np.asarray(g), 0, 'Sigmaz')
        J = np.asarray(J)
        if conserve is None:
            for u1, u2, dx in self.lat.nearest_neighbors:
                self.add_coupling(-J, u1, 'Sigmax', u2, 'Sigmax', dx)
        else:
            for op1, op2 in [('Sp', 'Sp'), ('Sp', 'Sm'), ('Sm', 'Sp'), ('Sm', 'Sm')]:
                for u1, u2, dx in self.lat.nearest_neighbors:
                    self.add_coupling(-J, u1, op1, u2, op2, dx)
        # 7) initialize MPO
        MPOModel.__init__(self, lat, self.calc_H_MPO())
        # skip 8): not a NearestNeighborModel...

In that way we can avoid generating multiple models for each of the classes, while keeping the current lattice interface including the sites.
We would only need to adjust the lattices classes Honeycomb and Kagome to accep a single "site" for the unit cell if it is the same for both/all three sites. We should probably do the latter anyways

[Johannes]

I followed your suggestions and moved the boundary conditions completely into the lattice, thanks for bringing that up!
As a side remark: I still distinguish the "lattice" boundary conditions bc strictly from the bc_MPS. This is necessary since somebody might want to consider a finite MPS on a chain, but use PBC for the couplings. Although this is not ideal for MPS, it can still be usefull for topoligical phases to enforce a unique groundstate...

The TFIsing2D is now indeed implemented as suggested above, allowing to select any of the predefined lattices.

I've also rewritten the Introduction to models, explaining a couple more things. If you (or anybody else) feels like there's still something unclear, let me now

[Johannes]

In response to the pull request #10 and to continue with this thread:

I completely agree that it's stupid to define each model multiple times. Using the trick above I got it down to two classes - one for a Chain derived from the NearestNeigborModel and once for generic 2D lattices (which include quasi-1D cases). Yet, I still define the couplings twice, which seems unnecessary and will lead to errors (e.g. when someone updates only one of the two classes).
The two classes are necessary as one is derived from the NearestNeigborModel, but defining the couplings twice isn't.
We should first define the model, allowing for general lattices (including a Chain), and then define the 1D version by deriving from it and the NearestNeigborModel.

I still like the idea to allow for a conserve='best' paramter, which checks the other model parameters on what can be preserved, as it done e.g. in the SpinChain model. This requires defining the Site instance inside the model, though, which makes everything a bit tricky in terms of where we initialize what. Personally, I also find it very convenient to have a single dictionary of parameters defining the model, including the lattice and what's conserved.

What do you think of changing the classes in tf_ising.py to the following:

Python: Select all

class TFIModel(CouplingModel, MPOModel):
    r"""Transverse field Ising model on a general (2D) lattice.
    ... documentation ..."""
    def __init__(self, model_param):
        # 0) read out/set default parameters
        lat = get_parameter(model_param, 'lattice', "Square", self.__class__)
        if isinstance(lat, str):
            Lx = get_parameter(model_param, 'Lx', 1, self.__class__)
            Ly = get_parameter(model_param, 'Ly', 4, self.__class__)
            bc_MPS = get_parameter(model_param, 'bc_MPS', 'infinite', self.__class__)
            bc_y = get_parameter(model_param, 'bc_y', 'cylinder', self.__class__)
            order = get_parameter(model_param, 'order', 'default', self.__class__)
            conserve = get_parameter(model_param, 'conserve', None, self.__class__)
            assert conserve != 'Sz'  # invalid!
            assert bc_y in ['cylinder', 'ladder']
            # 1-3)
            site = SpinHalfSite(conserve=conserve)
            # 4) lattice
            bc_x = 'periodic' if bc_MPS == 'infinite' else 'open'
            bc_y = 'periodic' if bc_y == 'cylinder' else 'open'
            lat = lattice.__dict__.get(lat)  # the corresponding class
            lat = lat(Lx, Ly, site, order=order, bc=[bc_x, bc_y], bc_MPS=bc_MPS)  # instance
        J = get_parameter(model_param, 'J', 1., self.__class__)
        g = get_parameter(model_param, 'g', 1., self.__class__)
        unused_parameters(model_param, self.__class__)  # checks for mistyped parameters
        # 5) initialize CouplingModel
        CouplingModel.__init__(self, lat)
        # 6) add terms of the Hamiltonian
        # (u is always 0 as we have only one site in the unit cell)
        self.add_onsite(-np.asarray(g), 0, 'Sigmaz')
        J = np.asarray(J)
        if all([site.conserve is None for site in lat.unit_cell]):
            for u1, u2, dx in self.lat.nearest_neighbors:
                self.add_coupling(-J, u1, 'Sigmax', u2, 'Sigmax', dx)
        else:
            for op1, op2 in [('Sp', 'Sp'), ('Sp', 'Sm'), ('Sm', 'Sp'), ('Sm', 'Sm')]:
                for u1, u2, dx in self.lat.nearest_neighbors:
                    self.add_coupling(-J, u1, op1, u2, op2, dx)
        # 7) initialize MPO
        MPOModel.__init__(self, lat, self.calc_H_MPO())


class TFIChain(TFIModel, NearestNeighborModel):
    r"""Transverse field Ising model on a chain. 
    ... documentation ... """
    def __init__(self, model_param):
        # 0) read out/set default parameters
        L = get_parameter(model_param, 'L', 2, self.__class__)
        bc_MPS = get_parameter(model_param, 'bc_MPS', 'finite', self.__class__)
        conserve = get_parameter(model_param, 'conserve', 'parity', self.__class__)
        assert conserve != 'Sz'
        # 1-3)
        site = SpinHalfSite(conserve=conserve)
        # 4) lattice
        bc = 'periodic' if bc_MPS == 'infinite' else 'open'
        lat = Chain(L, site, bc=bc, bc_MPS=bc_MPS)
        model_param['lattice'] = lat
        super().__init__(model_param)
        # 8) initialize bonds
        NearestNeighborModel.__init__(self, lat, self.calc_H_bond())

Note that this even keeps backwards compatibility
The general TFIModel is 2D in the sense that generating the lattice with lat(Lx, Ly, site, order=order, bc=[bc_x, bc_y], bc_MPS=bc_MPS) is the common call structure for 2D lattices, which doesn't work for a Chain, but it can also be used with any other, user-defined lattice if an instance is directly given with the lattice parameter.

All other "...Chain" models can be generalized in a similar way to 2D lattices, most importantly the SpinChain. I guess we should do that when we include this change into the master branch....

[Johannes]

I've updated the main repository with quite some changes in the models, so let me comment on this once more.

Making the jungle of classes even worse, I've added two more classes: a common base Model and the CouplingMPOModel. Moreover, I've update most of the models in TeNPy to use this new CouplingMPOModel.
The CouplingMPOModel follows exactly the idea outlined in the previous post, and suggested by [mention]mrader[/mention]:
We have a model parameter lattice. It can eather be directly an Instance of a lattice, or it can be a string like "Chain" (default value) or "Honeycomb". In the latter cases, further model parameters are read out (length in x,y direction, boundary conditions, ...) to initialize the corresponding lattice. In that way, the models which are defined as by (next-) nearest neighbor terms on an arbitrary lattice get all the same class.
However, some models need very specific lattices (e.g. the Kitaev-Honeycomb model is only defined on a Honeycomb lattice), in which case the init_model method might just be overwritten.

I've expanded the Introduction to models by another chapter, explaining how to write models with the CouplingMPOModel.

Another thing I've implemented is that one can now group neighboring sites (in the MPS sense). This is an example why we actually need the different classes for the different representations: depending on which representation(s) we have, we have to do different things for the grouping.
With the approach of the different classes, we can solve problem this very elegantly by overwriting the method group_sites and a clever use of super(). The common Model base class became necessary for this, but makes also sense regarding the common structure of defining self.lat.
Another example would be a function to