"""
    hubbardMPO(t, U, pchs, N[; T = ComplexF64])
returns a `U1`-symmetric MPO for the Hubbard model
```
    H = -∑_i t (a^†_i a_(i+1) + a_i a^†_(i+1)) + U/2 n_i (n_i-1)
```
on N sites.
"""
function hubbardMPO(t, U, pchs, N; T = ComplexF64)
    pdims = fill(1, length(pchs))
    bulk = DASTensor{T,4}(
        U1(),
        (U1Charges(-1:1), pchs, pchs, U1Charges(-1:1)),
        ([1,2,1], copy(pdims), copy(pdims), [1,2,1]),
        InOut(1,1,-1,-1))
    initwithzero!(bulk)
    preapp  = T[0 0; 0 1] # |0><0|_p ⊗ |1><1|_d
    postapp = T[1 0; 0 0] # |0><0|_p ⊗ |0><0|_d
    app     = T[0 0; 1 0] # |0><0|_p ⊗ |1><0|_d
    mat00to1 = reshape(T[0, 1], 2, 1)
    mat1to01 = reshape(T[1, 0], 1, 2)

    for u1ch in pchs
        ch = u1ch.ch
        peye= Matrix{T}(I, 1, 1)
        #local interaction term
        @tensor deg[1, 2, 3, 4] := preapp[1, 4]  * peye[2, 3] +  # |1><1| * Ι
                                   postapp[1, 4] * peye[2, 3] +  # |0><0| * Ι
                                   U/2 * ch * (ch-1)  * app[1, 4] * peye[2, 3]
                                   # U * ch * (ch-1) * |1><0|
        bulk[DASSector(U1Charge(0), U1Charge(ch), U1Charge(ch), U1Charge(0))] = deg

        #hopping term
        U1Charge(ch-1) in pchs || continue
        #a a^†
        # (0,0) → a_i → (1,0)
        @tensor deg[1, 2, 3, 4] := mat00to1[1, 4] * peye[2, 3]
        bulk[DASSector(U1Charge(0), U1Charge(ch), U1Charge(ch-1), U1Charge(1))] = -t*deg
        # (1,0) → a^†_{i+1} → (0,1)
        @tensor deg[1, 2, 3, 4] := mat1to01[1, 4] * peye[2, 3]
        bulk[DASSector(U1Charge(1), U1Charge(ch-1), U1Charge(ch), U1Charge(0))] = deg

        #a^† a
        # (0,0) → a^†_i → (-1,0)
        @tensor deg[1, 2, 3, 4] := mat00to1[1, 4] * peye[2, 3]
        bulk[DASSector(U1Charge(0), U1Charge(ch-1), U1Charge(ch), U1Charge(-1))] = -t*deg
        # (-1,0) → a_{i+1} → (0,1)
        @tensor deg[1, 2, 3, 4] := mat1to01[1, 4] * peye[2, 3]
        bulk[DASSector(U1Charge(-1), U1Charge(ch), U1Charge(ch-1), U1Charge(0))] = deg
    end
    # (3)        (4)       (3)
    # [lb](2) (2)[b](3) (2)[rb]
    # (1)        (1)       (1)
    lb = DASTensor{T,1}(U1(), (U1Charges(-1:1),), ([1,2,1],), InOut(-1))
    rb = DASTensor{T,1}(U1(), (U1Charges(-1:1),), ([1,2,1],), InOut( 1))
    lb[DASSector(U1Charge(0))] = [0, 1]
    rb[DASSector(U1Charge(0))] = [1, 0]

    @tensor bulk[1,2,3,4]       := bulk[2,1,4,3]
    @tensor leftbound[1, 2, 3]  := lb[-1] *  bulk[1, -1, 2, 3]
    @tensor rightbound[1, 2, 3] := bulk[1, 2, -1, 3] * rb[-1]

    return MPO{N}(leftbound, bulk,  rightbound)
end

"""
    creationOp(pchs::U1Charges [, T = ComplexF64])
return the creation operator `a†` acting on charges `pchs`
"""
function creationOp(pchs::U1Charges, T = ComplexF64)
    a = DASTensor{T,2}(U1(), (pchs, pchs),
        (fill(1, length(pchs)), fill(1, length(pchs))),
        InOut(1,-1))
    initwithzero!(a)
    for ch in pchs
        nch = ch ⊕ U1Charge(-1)
        nch in pchs || continue
        a[DASSector(nch, ch)] = reshape([sqrt(ch.ch)],1,1)
    end
    return a
end

"""
    annihilationOp(pchs::U1Charges [, T = ComplexF64])
return the annihilation operator `a` acting on charges `pchs`
"""
function annihilationOp(pchs::U1Charges, T = ComplexF64)
    a = DASTensor{T,2}(U1(), (pchs, pchs),
        (fill(1, length(pchs)), fill(1, length(pchs))),
        InOut(1,-1))
    initwithzero!(a)
    for ch in pchs
        nch = ch ⊕ U1Charge(-1)
        nch in pchs || continue
        a[DASSector(ch, nch)] = reshape([sqrt(ch.ch)],1,1)
    end
    return a
end
