\[ H = -t \sum_{< i,j >} \left( a_i^+ a_j + b_i^+ b_j + H.c. \right) + \frac{U}{2} \sum_i [a_i^+ a_i (a_i^+ a_i - 1) + b_i^+ b_i (b_i^+ b_i - 1)] + U_{AB} \sum_i a_i^+ a_i b_i^+ b_i, \]
where \(U\) and \(U_{AB}\) are the inter- and intrastate onsite interactions, respectively. I have implemented this as its own `CouplingMPOModel` class, while heeding some of the things I've learned from studying this topic on combining multiple charges on the same site.
Essentially I set up a superlattice with two sites in every unit cell: an _A_ site and a _B_ site. In the `init_lattice()` routine I use a `pairs` dictionary to get the geometry right:
Code: Select all
class TwoComponentBoseHubbardModel(CouplingMPOModel):
    def init_sites(self, model_params):
        n_max = model_params.get('n_max', 3)
        filling = model_params.get('filling', 1)
        conserve = model_params.get('conserve', 'N')
        if conserve == 'best':
            conserve = 'N'
            if self.verbose >= 1.:
                print(self.name + ": set conserve to", conserve)
        siteA = BosonSite(Nmax=n_max, conserve=conserve, filling=filling)
        siteB = BosonSite(Nmax=n_max, conserve=conserve, filling=filling)
        multi_sites_combine_charges([siteA, siteB])
        return [siteA, siteB]
        def init_lattice(self, model_params):
		pairs = {'onsite': [(0, 1, np.array([0]))], 'nearest_neighbors': [(0, 0, np.array([1])), (1, 1, np.array([1]))]}
		L = model_params.get('L', 10)
		sites = self.init_sites(model_params)
		lat = lattice.Lattice([L], sites, pairs = pairs)
		return lat	
	    def init_terms(self, model_params):
		t = model_params.get('t', 1.)
		U = model_params.get('U', 10.)
		UAB = model_params.get('UAB', 0.)
		mu = model_params.get('mu', 0)
		for u in range(len(self.lat.unit_cell)):
			self.add_onsite(-mu - U / 2., u, 'N')
			self.add_onsite(U / 2., u, 'NN')
		for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
			self.add_coupling(-t, u1, 'Bd', u2, 'B', dx, plus_hc=True)
		for u1, u2, dx in self.lat.pairs['onsite']:
			self.add_coupling(UAB, u1, 'N', u2, 'N', dx)
It does seem to me, though, that it runs a little inefficiently, so perhaps my choice of geometry is suboptimal. I need a bond dimension of about 300 to get a converged result for 20 sites. I must say I've never tried to use MPS algorithms on the Bose-Hubbard model before, so it may just be my inexperience (in which case, too bad!) On the other hand I wonder if the MPS representation is somehow more efficient if I use a ladder geometry, with the _A_ and _B_ sites being on either side of every rung.