Problems with creating multi-site coupling on a chain.

How do I use this algorithm? What does that parameter do?
Post Reply
frenzie92
Posts: 4
Joined: 14 Nov 2019, 16:07

Problems with creating multi-site coupling on a chain.

Post by frenzie92 »

Hello everyone.

So I want to implement the following spinless fermion Hamiltonian on a simple 1-D chain :

\(
H=\sum_i \left[ V_{10} n_i n_{i+1}+V_{20} n_i n_{i+2}-V_{21} \left( c^\dagger_i c_{i+1} c_{i+2} c^\dagger_{i+3} + h.c. \right)\right].
\)

In order to do this, I use add_multi_coupling as is done in the toric code example, and also mentioned in this topic.


What I've done so far is this (basically changed parts of the toric code example to fit my purposes. What I have is this

Code: Select all

import numpy as np

from .lattice import Lattice, Site, Chain
from ..networks.site import FermionSite
from .model import MultiCouplingModel, CouplingMPOModel
from ..tools.params import get_parameter
from ..tools.misc import any_nonzero

__all__ = ['FQHChain']

class FQHChain(CouplingMPOModel, MultiCouplingModel):
    r""" 

    Parameters
    ----------
    L : int
        Length of the chain.
    V10, V20, V21 : float
        Couplings as defined for the Hamiltonian above.
    bc_MPS : {'finite' | 'infinte'}
        MPS boundary conditions. Coupling boundary conditions are chosen appropriately.
    """

    def __init__(self, model_params):
        CouplingMPOModel.__init__(self, model_params)

    def init_sites(self, model_params):
        site = FermionSite(conserve='N')
        return site

    def init_lattice(self, model_params):
        site = self.init_sites(model_params)
        L = get_parameter(model_params, 'L', 2, self.name) # default values of parameters are currently arbitrary.
        bc_MPS = get_parameter(model_params, 'bc_MPS', 'finite', self.name)
        bc = 'periodic' if bc_MPS == 'finite' else 'open'
        lat = Chain(L, site, bc=bc, bc_MPS=bc_MPS) 
        return lat

    def init_terms(self, model_params):
        V10 = get_parameter(model_params, 'V10', 1., self.name, asarray=True)  # default values of parameters are currently arbitrary. 
        V20 = get_parameter(model_params, 'V20', 0.2, self.name, True)
        V21 = get_parameter(model_params, 'V21', 0.5, self.name, True)
        self.add_coupling(V10, 0, 'N', 0, 'N', 1)
        self.add_coupling(V20, 0, 'N', 0, 'N', 2)
        self.add_multi_coupling(-V21, 0, 'Cd', [(0, 'C', 1), (0, 'C', 2), (0, 'Cd', 3)] )
        self.add_multi_coupling(-V21, 0, 'C', [(0, 'Cd', -1), (0, 'Cd', -2), (0, 'C', -3)] )  # h.c. of previous term
        # done
I call this fqh_chain.py and save it in the same directory as the other models.

In order to check it, I run a variant of the code found in page 20 of the tenpy paper. In particular, I run the following

Code: Select all

import numpy as np
import random

from tenpy.networks.mps import MPS
from tenpy.models.fqh_chain import FQHChain
from tenpy.algorithms import dmrg


L=16
N=8

places=random.sample(range(1,L), N)


#We create a random state with N electrons in occupation number representation

initial_psi=['empty']*L


for i in range(0,len(places)):
	initial_psi[places[i]]='full'

M = FQHChain({"L": L, "V10": 1, "V20": 0.1, "V21": 0., "bc_MPS": "finite"})
psi = MPS.from_product_state(M.lat.mps_sites(), initial_psi, "finite")
dmrg_params = {"trunc_params": {"chi_max": 30, "svd_min": 1.e-10}}
dmrg.run(psi, M, dmrg_params) # find the ground state
print("E =", sum(psi.expectation_value(M.H_bond[1:])))
print("final bond dimensions: ", psi.chi)

The output is this

Code: Select all

parameter 'L'=16 for FQHChain
parameter 'bc_MPS'='finite' for FQHChain
parameter 'V10'=1 for FQHChain
parameter 'V20'=0.1 for FQHChain
parameter 'V21'=0.0 for FQHChain
parameter 'trunc_params'={'chi_max': 30, 'svd_min': 1e-10} for Sweep
================================================================================
sweep 10, age = 16
Energy = 3.3000000000000003, S = 0.0000000000000000, norm_err = 0.0e+00
Current memory usage 47704.0 MB, time elapsed: 3.0 s
Delta E = nan, Delta S = nan (per sweep)
max_trunc_err = 0.0000e+00, max_E_trunc = 4.4409e-16
MPS bond dimensions: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
================================================================================
sweep 20, age = 16
Energy = 3.3000000000000003, S = 0.0000000000000000, norm_err = 0.0e+00
Current memory usage 47860.0 MB, time elapsed: 6.0 s
Delta E = 0.0000e+00, Delta S = 0.0000e+00 (per sweep)
max_trunc_err = 0.0000e+00, max_E_trunc = 4.4409e-16
MPS bond dimensions: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
================================================================================
DMRG finished after 20 sweeps.
total size = 16, maximum chi = 1
================================================================================
Traceback (most recent call last):
  File "test.py", line 27, in <module>
    print("E =", sum(psi.expectation_value(M.H_bond[1:])))
AttributeError: 'FQHChain' object has no attribute 'H_bond'
There are a number of things wrong with this. First and foremost, the error AttributeError: 'FQHChain' object has no attribute 'H_bond'. Also, something's fishy with the bond dimensions, they're all unity! \(\Delta E =0\), and in general it looks like my attempt at implementing my own model failed miserably, which makes me very sad :cry: :cry:

To try to find out what's wrong I tried adding a similar four-site coupling to the xxz chain (in particular, XXZChain2). The Hamiltonian becomes

\(
H = \sum_i \left[ \frac{\mathtt{J_{xx}}}{2} (S^{+}_i S^{-}_{i+1} + S^{-}_i S^{+}_{i+1})
+ \mathtt{J_z} S^z_i S^z_{i+1}
- \sum_i \mathtt{h_z} S^z_i\right]+K \sum_i \left(S^z_i S^z_{i+1} S^z_{i+2} S^z_{i+3} + h.c. \right)
\)

So, I just add with add_multi_coupling the extra terms. I have

Code: Select all

import numpy as np

from .lattice import Site, Chain
from .model import CouplingModel, NearestNeighborModel, MPOModel, CouplingMPOModel, MultiCouplingModel
from ..linalg import np_conserved as npc
from ..tools.params import get_parameter, unused_parameters
from ..networks.site import SpinHalfSite  # if you want to use the predefined site

__all__ = ['TestChain']


class TestChain(CouplingMPOModel, NearestNeighborModel, MultiCouplingModel):
    """Another implementation of the Spin-1/2 XXZ chain with Sz conservation.

    This implementation takes the same parameters as the :class:`XXZChain`, but is implemented
    based on the :class:`~tenpy.models.model.CouplingMPOModel`.
    """
    def __init__(self, model_params):
        model_params.setdefault('lattice', "Chain")
        CouplingMPOModel.__init__(self, model_params)

    def init_sites(self, model_params):
        return SpinHalfSite(conserve='Sz')  # use predefined Site

    def init_terms(self, model_params):
        # read out parameters
        Jxx = get_parameter(model_params, 'Jxx', 1., self.name, True)
        Jz = get_parameter(model_params, 'Jz', 1., self.name, True)
        hz = get_parameter(model_params, 'hz', 0., self.name, True)
        Jv = get_parameter(model_params, 'Jv', 0., self.name, True)
        # add terms
        self.add_onsite(-hz, 0, 'Sz')
        self.add_coupling(Jxx * 0.5, 0, 'Sp', 0, 'Sm', 1)
        self.add_coupling(np.conj(Jxx * 0.5), 0, 'Sp', 0, 'Sm', -1)  # h.c.
        self.add_coupling(Jz, 0, 'Sz', 0, 'Sz', 1)
        self.add_multi_coupling(Jv, 0, 'Sz', [(0, 'Sz', 1), (0, 'Sz', 2), #this is the only thing I add, along with the inclusion of MultiCouplingModel
                                                  (0, 'Sz', 3)])
        self.add_multi_coupling(Jv, 0, 'Sz', [(0, 'Sz', -1), (0, 'Sz', -2),
                                                  (0, 'Sz', -3)])
Let's call this 'test_chain.py'.

To test it run (note that I have set \(K=0\)).

Code: Select all

from tenpy.networks.mps import MPS
from tenpy.models.test_chain import TestChain
from tenpy.algorithms import dmrg


M = TestChain({"L": 4, "Jxx": 1, "Jz": 0.1, "hz": 0., "K":0, "bc_MPS": "finite"})
psi = MPS.from_product_state(M.lat.mps_sites(), [0,1,1,0], "finite")
dmrg_params = {"trunc_params": {"chi_max": 30, "svd_min": 1.e-10}}
dmrg.run(psi, M, dmrg_params) # find the ground state
print("E =", sum(psi.expectation_value(M.H_bond[1:])))
print("final bond dimensions: ", psi.chi)
This works quite well, giving

Code: Select all

parameter 'lattice'='Chain' for TestChain
parameter 'bc_MPS'='finite' for TestChain
parameter 'L'=4 for TestChain
parameter 'Jxx'=1 for TestChain
parameter 'Jz'=0.1 for TestChain
parameter 'hz'=0.0 for TestChain
parameter 'K'=0 for TestChain
parameter 'trunc_params'={'chi_max': 30, 'svd_min': 1e-10} for Sweep
================================================================================
sweep 10, age = 4
Energy = -1.1636015623654488, S = 0.5927718553565784, norm_err = 5.1e-16
Current memory usage 47788.0 MB, time elapsed: 0.6 s
Delta E = nan, Delta S = nan (per sweep)
max_trunc_err = 0.0000e+00, max_E_trunc = 1.1102e-15
MPS bond dimensions: [2, 4, 2]
================================================================================
sweep 20, age = 4
Energy = -1.1636015623654492, S = 0.5927718553565789, norm_err = 5.4e-16
Current memory usage 47836.0 MB, time elapsed: 1.1 s
Delta E = -4.4409e-17, Delta S = 5.5511e-17 (per sweep)
max_trunc_err = 0.0000e+00, max_E_trunc = 2.2204e-16
MPS bond dimensions: [2, 4, 2]
================================================================================
DMRG finished after 20 sweeps.
total size = 4, maximum chi = 4
================================================================================
E = -1.163601562365449
final bond dimensions:  [2, 4, 2]
Now try giving \(K\) a finite value, say \(K=0.01\), I get

Code: Select all

parameter 'lattice'='Chain' for TestChain
parameter 'bc_MPS'='finite' for TestChain
parameter 'L'=4 for TestChain
parameter 'Jxx'=1 for TestChain
parameter 'Jz'=0.1 for TestChain
parameter 'hz'=0.0 for TestChain
parameter 'K'=0.01 for TestChain
Traceback (most recent call last):
  File "test.py", line 9, in <module>
    M = TestChain({"L": 4, "Jxx": 1, "Jz": 0.1, "hz": 0., "K":0.01, "bc_MPS": "finite"})
  File "/Users/USERNAME/TeNPy/tenpy/models/test_chain.py", line 27, in __init__
    CouplingMPOModel.__init__(self, model_params)
  File "/Users/USERNAME/TeNPy/tenpy/models/model.py", line 1240, in __init__
    NearestNeighborModel.__init__(self, lat, self.calc_H_bond())
  File "/Users/USERNAME/TeNPy/tenpy/models/model.py", line 918, in calc_H_bond
    H_bond = ct.to_nn_bond_Arrays(sites)
  File "/Users/USERNAME/TeNPy/tenpy/networks/terms.py", line 635, in to_nn_bond_Arrays
    H_add = strength * npc.outer(site_i.get_op(op1), site_j.get_op(op2))
  File "/Users/USERNAME/TeNPy/tenpy/networks/site.py", line 305, in get_op
    names = name.split(' ')
AttributeError: 'tuple' object has no attribute 'split'
So now I get a different kind of error :? I noticed that not including the NearestNeighborModel yields the same 'AttributeError: 'TestChain' object has no attribute 'H_bond''. So I guess that's what the problem was initially. But, looking at the second model now, I still cannot have multi-site couplings without yielding errors :?: In particular, I get 'AttributeError: 'tuple' object has no attribute 'split''.

Sorry if this question was too long, I wanted to provide as much information as I could, and provide everything I've tried. So the basic problem, how do I get the initial model to give sensible results? What's wrong with the way I did it? Also, what's wrong with the second model (perhaps answering this might help with the first model - since the only difference with the already made XXZChain2 model is the inclusion of add_multi_coupling. If this is done I can rewrite my own model based on XXZChain2. )

Any help would be very much appreciated.
Last edited by frenzie92 on 21 Dec 2019, 15:35, edited 1 time in total.
User avatar
Johannes
Site Admin
Posts: 413
Joined: 21 Jul 2018, 12:52
Location: TU Munich

Re: Problems with creating multi-site coupling on a chain.

Post by Johannes »

By definition, a tenpy.models.model.NearestNeighborModel can only contain nearest-neighbor couplings.
That is the case because it provides a representation of the hamiltonian in terms of H_bond, where each H_bond is a term involving just two neighboring sites.
Hence, you can *not* make a model inherited from the NearestNeighborModel and add couplings over more than two neighboring sites. (You also cannot add two-site couplings with dx > 1.)

The correct way to go is hence to *not* inherit from the NearestNeighborModel, as you do in your FQHChain.
However, that implies that there also is no H_bond representation of the Hamiltonian, which you try to use in evaluating the energy.
It is an tenpy.models.model.MPOModel, however (inherited from the CouplingMPOModel[/tenpy], which means that it has a representation of the hamiltonian as an MPO, H_MPO.
Try to evaluate the energy with M.H_MPO.expectation_value(psi) instead, or even simipler just use the "E" returned by the DMRG run.

That should answer why the error appears in your first example, and why the second example fails.
I agree that the error message in the second case is quite obscure, though, and does not really tell you about the problem.
That was indeed a bug, and I fixed it now (in the newest version only, not yet in v0.5.0) to give a more sensible error message.


To your initial FQHChain model: that's fine!
The reason why you don't get anything else but a product state from DMRG is simply that you conserve the total particle number and initialize it in the vacuum state, i.e in the charge sector with 0 particles.
This charge sector consists only of that single vacuum state. DMRG is not allowed to change the charge sector, so it will stay there, and correctly output the vacuum state (which also is an Eigenstate of your hamiltonian, with eigenenergy 0).
The way to go is to choose a different initial state, e.g.

Code: Select all

initial_psi = ['empty', 'full']*(L//2)  # for even L
Also, you should put mixer=True in your DMRG parameters, to allow getting out of this initial state.
Post Reply