Particle Number Conservation

How do I use this algorithm? What does that parameter do?
Post Reply
QS Wang
Posts: 8
Joined: 22 Mar 2021, 10:27

Particle Number Conservation

Post by QS Wang »

Hello,

I am trying to find the ground state and calculate the ground state energy of a spinless fermi system using GroupedSite and DMRG. However, the problem is that I need the number of impurity to be fixed, while actually it always change and finally DMRG always give me the global minimum with another number of impurity and number of host particles. For example, when I set the number of host fermions and of the impurity to be 10 and 1 respectively in the initial state, DMRG would give me the ground state with 5.5 fermions and 5.5 impurities, no matter the impurity is fermion or boson. How could I keep the number of host particles and impurities conserved during DMRG?

Thanks.

The definition of the model is

Code: Select all

#define the model with grouped sites for fermionic impurity

class GroupedBosonImpModel(CouplingMPOModel):

  def init_sites(self, model_params):
    conserve_F = model_params.get('conserve_F', 'N')
    filling_F = model_params.get('filling_F',0.5)
    conserve_B = model_params.get('conserve_B', 'N')
    filling_B = model_params.get('filling_B',0)
    Nmax = model_params.get('Nmax',1)
    model_label = model_params.get('label',None)
    model_charges = model_params.get('charges','same')

    HostSite = FermionSite(conserve_F,filling_F)
    ImpSite = BosonSite(Nmax,conserve_B,filling_B)

    site = GroupedSite([HostSite,ImpSite])
    NH=site.N0
    NI=site.N1
    N_product=npc.tensordot(NH,NI,[1,0])
    site.add_op('N_product',N_product)
    return site
  
  
  def init_terms(self, model_params):

    #confirm that the input lattice is a chain
    if not isinstance(self.lat,Chain):
      raise ValueError("The input lattice is not a Chain") 

    #get the coupling constants, t*c_{i+1}^\dagger c_i, 
    # h*x_{i+1}^\dagger x_i
    # V*n_e*n_e and g*n_x*n_e
    # mue, mux the potential for particle number conservation

    t = model_params.get('t', 1.)
    h = model_params.get('h', 1.)
    V = model_params.get('V', 0.)
    g = model_params.get('g', 0.)
    mue = model_params.get('mue', 1.)
    mux = model_params.get('mux', 1.)

    #implement chemical potential
    self.add_onsite(-mue,0,'N0',plus_hc=True) 
    self.add_onsite(-mux,0,'N1',plus_hc=True)

    #implement e-imp interaction
    self.add_onsite(-g,0,'N_product',plus_hc=True) 

    #implement e-e interaction
    self.add_coupling(-V,0,'dN0',0,'dN0',1,plus_hc=True)

    #implement kinetic energy
    self.add_coupling(-t,0,'Cd0',0,'C0',1,plus_hc=True)
    self.add_coupling(-h,0,'Bd1',0,'B1',1,plus_hc=True)

The code for dmrg is

Code: Select all

p_state=[[2],[0]]
for i in range(1,N):
  p_state=p_state+[[2],[0]]

import time
start_time = time.time()

psi_trial=MPS.from_lat_product_state(TestModel.lat,p_state)
Nimp_0=psi_trial.expectation_value('N1')
Ne_0=psi_trial.expectation_value('N0')
print('Nimp_0=',Nimp_0.sum())
print('Ne_0=',Ne_0.sum())

dmrg_params_groupedsite = {
    'mixer': None,  # setting this to True helps to escape local minima
    'max_E_err': 1.e-10,
    'trunc_params': {
        'chi_max': 10,
        'svd_min': 1.e-10,
    },
    'verbose': True,
    'combine': True
}

dmrg_params = {
    'mixer': None,
    'max_E_err': 1.e-10,
}
eng_GS = dmrg.TwoSiteDMRGEngine(psi_trial, TestModel, dmrg_params_groupedsite)
E, psi = eng_GS.run() # the main work; modifies psi in place

print("--- %s seconds ---" % (time.time() - start_time))

print('E=', E)


Nimp=psi.expectation_value('N1')
Ne=psi.expectation_value('N0')

print('Nimp=',Nimp.sum())
print('Ne=',Ne.sum())
and the output is

Nimp_0= 0.0
Ne_0= 10.0
final DMRG state not in canonical form up to norm_tol=1.00e-05: norm_err=2.08e-02
--- 10.078522682189941 seconds ---
E= -40.299434837756976
Nimp= 5.000000000000001
Ne= 5.000000000000001
User avatar
Johannes
Site Admin
Posts: 413
Joined: 21 Jul 2018, 12:52
Location: TU Munich

Re: Particle Number Conservation

Post by Johannes »

Your couplings preserve the number of fermions and bosons independently, so you want to end up with qnumber=2 charges. Try
using charges='independent' when combining the sites.

Moreover, you double-count the N-terms compared to the hopping terms by explicitly adding plus_hc=True, at least if you intended to implement the following Hamiltonian
\[ H = - t \sum_j ((c^e_j)^\dagger c^e_{j+1} + h.c.) - h \sum_j ((b^i_j)^\dagger b^i_{j+1} + h.c.) - V \sum_j dN^e_j dN^e_{j+1} - g \sum_j N^e_j N^i_j - \mu^e \sum_j N^e_j - \mu^i \sum_j N^i_j \]

That being said, there is actually no reason to combine the charges here, and doing so will make DMRG run slower by quite a bit when you go to larger system sizes, in particular when your N_max is not just 2.
Try the following instead:

Code: Select all

from tenpy.networks.site import set_common_charges
import tenpy

#define the model with grouped sites for fermionic impurity
class GroupedBosonImpModel(CouplingMPOModel):
  default_lattice = "Ladder"
  force_default_lattice = True
  # this requires TeNPy v0.8.0

  def init_sites(self, model_params):
    conserve_F = model_params.get('conserve_F', 'N')
    filling_F = model_params.get('filling_F',0.5)
    conserve_B = model_params.get('conserve_B', 'N')
    filling_B = model_params.get('filling_B',0)
    Nmax = model_params.get('Nmax',1)
    model_label = model_params.get('label',None)
    model_charges = model_params.get('charges','independent')

    HostSite = FermionSite(conserve_F,filling_F)
    ImpSite = BosonSite(Nmax,conserve_B,filling_B)
    set_common_charges([HostSite, ImpSite], model_charges)
    return [HostSite, ImpSite]
  
  def init_terms(self, model_params):

    #confirm that the input lattice is a chain
    if not isinstance(self.lat, tenpy.models.lattice.Ladder):
      raise ValueError("The input lattice is not a Ladder") 

    #get the coupling constants, t*c_{i+1}^\dagger c_i, 
    # h*x_{i+1}^\dagger x_i
    # V*n_e*n_e and g*n_x*n_e
    # mue, mux the potential for particle number conservation

    t = model_params.get('t', 1.)
    h = model_params.get('h', 1.)
    V = model_params.get('V', 0.)
    g = model_params.get('g', 0.)
    mue = model_params.get('mue', 1.)
    mux = model_params.get('mux', 1.)

    #implement chemical potential
    self.add_onsite(-mue,0,'N')  # `u=0` = HostSite
    self.add_onsite(-mux,1,'N')  # `u=1` = ImpSite

    #implement e-imp interaction
    self.add_coupling(-g, 0, 'N', 1, 'N', dx=[0])

    #implement e-e interaction
    self.add_coupling(-V,0,'dN',0,'dN',dx=[1])

    #implement kinetic energy
    self.add_coupling(-t,0,'Cd0',0,'C0',1,plus_hc=True)
    self.add_coupling(-h,0,'Bd1',0,'B1',1,plus_hc=True)
Your inintial state would then be something like

Code: Select all

Lx = M.lat.Ls[0]
assert Lx % 2 == 0
p_state = [['full', '1'], ['empty', '0']] * (Lx//2)
QS Wang
Posts: 8
Joined: 22 Mar 2021, 10:27

Re: Particle Number Conservation

Post by QS Wang »

All of plus_hc are set to be True since tenpy always raises TypeError once it is False.

Code: Select all

TypeError                                 Traceback (most recent call last)
<ipython-input-5-73d291e71b76> in <module>()
     18 }
     19 
---> 20 TestModel=GroupedBosonImpModel(test_params)

3 frames
/usr/local/lib/python3.7/dist-packages/tenpy/models/model.py in __init__(self, model_params)
   1711         CouplingModel.__init__(self, lat, explicit_plus_hc=explicit_plus_hc)
   1712         # 6) add terms of the Hamiltonian
-> 1713         self.init_terms(model_params)
   1714         # 7-8) initialize H_MPO, and H_bonds, if necessary
   1715         self.init_H_from_terms()

<ipython-input-3-dd0d33b907cb> in init_terms(self, model_params)
     45 
     46     #implement e-imp interaction
---> 47     self.add_coupling(-g, 0, 'N', 1, 'N', dx=[0])
     48 
     49     #implement e-e interaction

/usr/local/lib/python3.7/dist-packages/tenpy/models/model.py in add_term_function(self, *args, **kwargs)
   1624     @wraps(f)
   1625     def add_term_function(self, *args, **kwargs):
-> 1626         res = f(self, *args, **kwargs)
   1627         if hasattr(self, 'H_MPO') and not getattr(self, 'manually_call_init_H', False):
   1628             warnings.warn(

/usr/local/lib/python3.7/dist-packages/tenpy/models/model.py in add_coupling(self, strength, u1, op1, u2, op2, dx, op_string, str_on_first, raise_op2_left, category, plus_hc)
   1007                 # (this reduces the MPO bond dimension with `explicit_plus_hc=True`)
   1008             else:
-> 1009                 strength_vals /= 2.  # ... so we should avoid double-counting
   1010         if category is None:
   1011             category = "{op1}_i {op2}_j".format(op1=op1, op2=op2)

TypeError: ufunc 'true_divide' output (typecode 'd') could not be coerced to provided output parameter (typecode 'l') according to the casting rule ''same_kind''
I have been confused with this problem for days and finally decided to correct the coefficients after the simulation by hand.

Thanks a lot.
User avatar
Johannes
Site Admin
Posts: 413
Joined: 21 Jul 2018, 12:52
Location: TU Munich

Re: Particle Number Conservation

Post by Johannes »

As far as I understand the error message, it claims that strength_vals is an integer array and hence can not divided by 2. in place. On your side, there's a very easy fix: specifying g = 1. as a float instead of g = 1 as an int - as you do in the default parameters.
On TeNPy side, this is likely something people will run into again, so I just fixed it in 050d311. Unfortunately, it didn't make it into version 0.8.0 any more, which I released last Friday...
QS Wang
Posts: 8
Joined: 22 Mar 2021, 10:27

Re: Particle Number Conservation

Post by QS Wang »

After attaining the ground state I tried to apply \(C^\dagger_p\) on it to calculate the energy spectrum. To do so I need to apply Cd to every sites and times the result with some phase factor. However, I was met with a error saying "s has wrong shape: (1,) instead of 0" (shown below).

The code is

Code: Select all

psi0=psi
psi0.apply_local_op(1,'N')
and the result is

Code: Select all

ValueError                                Traceback (most recent call last)
<ipython-input-52-119ce83ffbf1> in <module>()
      1 psi0=psi
----> 2 psi0.apply_local_op(1,'N')

6 frames
/usr/local/lib/python3.7/dist-packages/tenpy/networks/mps.py in apply_local_op(self, i, op, unitary, renormalize, cutoff)
   2919                 self.set_SR(i + j, split_th._S[j + 1])
   2920         if not unitary:
-> 2921             self.canonical_form(renormalize)
   2922 
   2923     def apply_product_op(self, ops, unitary=None, renormalize=False):

/usr/local/lib/python3.7/dist-packages/tenpy/networks/mps.py in canonical_form(self, renormalize)
   2551         """
   2552         if self.finite:
-> 2553             return self.canonical_form_finite(renormalize)
   2554         else:
   2555             self.canonical_form_infinite(renormalize)

/usr/local/lib/python3.7/dist-packages/tenpy/networks/mps.py in canonical_form_finite(self, renormalize, cutoff)
   2602         else:
   2603             # we actually had a canonical form before, so we should *not* ignore the 'S'
-> 2604             M = self.get_B(0, form='Th')
   2605             form = 'B'  # for other 'M'
   2606         if self.bc == 'segment':

/usr/local/lib/python3.7/dist-packages/tenpy/networks/mps.py in get_B(self, i, form, copy, cutoff, label_p)
    878                 B = self._scale_axis_B(B, self.get_SL(i), new_form[0] - old_form[0], 'vL', cutoff)
    879             if new_form[1] is not None and new_form[1] - old_form[1] != 0.:
--> 880                 B = self._scale_axis_B(B, self.get_SR(i), new_form[1] - old_form[1], 'vR', cutoff)
    881         if label_p is not None:
    882             B = self._replace_p_label(B, label_p)

/usr/local/lib/python3.7/dist-packages/tenpy/networks/mps.py in _scale_axis_B(self, B, S, form_diff, axis_B, cutoff)
   3409             elif form_diff != 1.:
   3410                 S = S**form_diff
-> 3411             return B.scale_axis(S, axis_B)
   3412         else:
   3413             # e.g. during DMRG with a DensityMatrixMixer

/usr/local/lib/python3.7/dist-packages/tenpy/linalg/np_conserved.py in scale_axis(self, s, axis)
   2008         res = self.copy(deep=False)
   2009         res._qdata = res._qdata.copy()
-> 2010         res.iscale_axis(s, axis)
   2011         return res
   2012 

/usr/local/lib/python3.7/dist-packages/tenpy/linalg/np_conserved.py in iscale_axis(self, s, axis)
   1989         if s.shape != (self.shape[axis], ):
   1990             raise ValueError("s has wrong shape: " + str(s.shape) + " instead of " +
-> 1991                              str(self.shape[axis]))
   1992         self.dtype = np.find_common_type([self.dtype], [s.dtype])
   1993         leg = self.legs[axis]

ValueError: s has wrong shape: (1,) instead of 0
I tried it on groupedsite model as well and same error arise. It is confused as yesterday it still works quite well as I remember.

Sorry for the disturbance.
User avatar
Johannes
Site Admin
Posts: 413
Joined: 21 Jul 2018, 12:52
Location: TU Munich

Re: Particle Number Conservation

Post by Johannes »

It seems there's something gone wrong beforehand; if you try psi.test_sanity() before the psi0.apply_local_op(...), it will likely give you an error, right?
If this code is inside a loop, make sure you have actual copies psi0 = psi.copy(), not just new references.

That being said: in what way can you get the energy spectrum from applying C^\dagger_p on each site? What exactly are you trying to calculate?
QS Wang
Posts: 8
Joined: 22 Mar 2021, 10:27

Re: Particle Number Conservation

Post by QS Wang »

What I want to calculate is \(\langle GS|C_p H C_p^\dagger |GS \rangle \). I planed to apply Cd on each site and sum them up with some phase factors to construct \(C_p^\dagger GS \rangle \) and calculate the expectation of \(H\) for this state.
QS Wang
Posts: 8
Joined: 22 Mar 2021, 10:27

Re: Particle Number Conservation

Post by QS Wang »

Hi, it seems that this problem only appear at the first site (or the first two sites for ladder). That is, for example, if I am working with a ladder with length of 10, then the problem will arise at site 0 and up to now for i>1 the function MPS.apply_local_op(i,op) works well.
Post Reply