Non-Hermitian problem while using self-defined model

How do I use this algorithm? What does that parameter do?
Post Reply
mount
Posts: 2
Joined: 05 Mar 2021, 07:49

Non-Hermitian problem while using self-defined model

Post by mount »

Hi,

I am trying to define the Anyon-Hubbard model(AHM):

\(H = -J {\sum}^{L}_{j=2} ({\hat b}^{\dagger}_{j} {\hat b}_{j-1} e^{i \theta {\hat n}_{j}} + H.c.) + U \sum^{L}_{j=1} {\hat n}_{j} ({\hat n}_{j} - 1)\)

Code: Select all

import numpy as np
import matplotlib.pyplot as plt

import tenpy
import tenpy.linalg.np_conserved as npc
from tenpy.algorithms import dmrg
from tenpy.networks.mps import MPS

from tenpy.networks.site import BosonSite
from tenpy.models.lattice import Chain
from tenpy.models.model import CouplingModel, NearestNeighborModel, MPOModel

"""Definition of a model: the Anyon-Hubbard chain."""
class AHChain(CouplingModel, NearestNeighborModel, MPOModel):
    def __init__(self, L=20, Nmax=2, filling=0.2, J=1, U=1, theta=0.5*np.pi):
        boson = BosonSite(Nmax=Nmax, conserve='N', filling=filling)
        boson.add_op('gJWd', npc.expm(-1.j*theta*boson.N), hc='gJW')
        boson.add_op('gJW', npc.expm(1.j*theta*boson.N), hc='gJWd')
        # the lattice defines the geometry
        lattice = Chain(L, boson, bc='open', bc_MPS='finite')
        CouplingModel.__init__(self, lattice)
        # add terms of the Hamiltonian
        self.add_coupling(-J, 0, 'Bd', 0, 'B', -1, 'gJW', str_on_first=True, plus_hc=True)
        self.add_onsite(U, 0, 'NN')
        self.add_onsite(-U, 0, 'N')
        
        # finish initialization
        # generate MPO for DMRG
        MPOModel.__init__(self, lattice, self.calc_H_MPO())
        # generate H_bond for TEBD
        NearestNeighborModel.__init__(self, lattice, self.calc_H_bond())
However, it seems the Hamiltonian of the AHM is considered a non-Hermitian Hamiltonian:

Code: Select all

M = AHChain(L=20, Nmax=2, filling=0.2, J=1, U=0, theta=0.5*np.pi)
p_state = [[0], [0], [0], [0], [1]] * 4
psi = MPS.from_lat_product_state(M.lat, p_state)
The output of `M.H_MPO.is_hermitian()` is `False`

And if I try to run the DMRG algorithm:

Code: Select all

dmrg_params = {
    'mixer': None,  # setting this to True helps to escape local minima
    'max_E_err': 1.e-8,
    'trunc_params': {
        'chi_max': 100,
        'svd_min': 1.e-8,
    },
    'verbose': True,
    'combine': True
}
eng = dmrg.TwoSiteDMRGEngine(psi, M, dmrg_params) 
E, psi = eng.run() # the main work; modifies psi in place
Two warnings appear:

`...\Anaconda3\lib\site-packages\tenpy\algorithms\dmrg.py:2032: UserWarning: H is zero in the given block, nothing to diagonalize.We just return the initial state again.
warnings.warn("H is zero in the given block, nothing to diagonalize."
...\Anaconda3\lib\site-packages\tenpy\linalg\lanczos.py:250: UserWarning: Poorly conditioned Lanczos!
warnings.warn("Poorly conditioned Lanczos!")`

Finally, I get another UserWarning:

`...\Anaconda3\lib\site-packages\tenpy\algorithms\dmrg.py:420: UserWarning: final DMRG state not in canonical form within 'norm_tol' = 1.00e-05
warnings.warn(msg.format(nt=norm_tol))`

Of course, the result is not correct... Any help will be appreciated!

Thanks.
Last edited by mount on 12 Mar 2021, 02:01, edited 1 time in total.
User avatar
Johannes
Site Admin
Posts: 413
Joined: 21 Jul 2018, 12:52
Location: TU Munich

Re: Non-Hermitian problem while using self-defined model

Post by Johannes »

Wow, that sounds like some interesting physics!

Sometimes, it's helpful to double-check the terms produced by the model. In your case, print(M.all_coupling_terms().to_TermList()) prints

Code: Select all

-1.00000 * B gJW_0 Bd_1 +
-1.00000 * Bd gJWd_0 B_1 +
-1.00000 * B gJW_1 Bd_2 +
-1.00000 * Bd gJWd_1 B_2 +
...
Here, "B gJW_0" are the operators "B" and "gJW" in that order acting on site 0. The correct hermitian conjugate would be
"gJWd Bd", but that's not what you find in line 2. Now, you can ask: why?
Digging into the code of add_coupling, we find that there is a hard-coded minus sign in the case i > j.
Stricktly speaking, you might argue that this is a bug, so I just fixed it in f789cac8e365314310b304ef34abac1989b54ff3.
As a warning, let me note that there is other code where there is a hard-coded minus sign in case operators which require a JW string are exchanged, e.g. in tenpy.networks.terms.order_combine_term.
After that fixing commit, the above print now yields:

Code: Select all

-1.00000 * Bd gJWd_0 B_1 +
-1.00000 * gJW B_0 Bd_1 +
-1.00000 * Bd gJWd_1 B_2 +
-1.00000 * gJW B_1 Bd_2 +
...
That looks better :)

Nevertheless, I believe, your code will still not do what you really expect. As you might have noted from above, the gJW was added to the left site instead of the right (as you wrote it in your hamiltonian). Or did you actually mean to write `j-1` in the exponent?

Just for reference, I think one key information you might have missed is that you can simply add the operator names separated by a space to multiply them on-site. In other words, instead of using the op_str to add the term on one of them, you can also do this:

Code: Select all

self.add_coupling(-J, 0, 'Bd gJW', 0, 'B', -1, plus_hc=True)
However, using the "op_str" might be "more" correct in case you want to generalize this to longer-range hopping.
A clearer solution might be to actually define your own "AnionSite" with a non-trivial `JW` operator.
You might run into issues with `JW` vs `dagger(JW)` in several places, that would need to be fixed, though - I didn't think of such a non-trivial JW when I wrote the code.
mount
Posts: 2
Joined: 05 Mar 2021, 07:49

Re: Non-Hermitian problem while using self-defined model

Post by mount »

Hi Johannes,

Thanks a lot for your help!

I try this code:

Code: Select all

self.add_coupling(-J, 0, 'Bd gJW', 0, 'B', -1, plus_hc=True)
and it works well.

Actually, I wanted to add \(\hat n_j\) in the exponent at first instead of \(\hat n_{j-1}\). However, what I would like to mention is that \(\hat n_j\) and \(\hat n_{j-1}\) give the same physics...haha
Post Reply