from quspin.operators import hamiltonian # Hamiltonians and operators
from quspin.basis import spinless_fermion_basis_1d, tensor_basis # Hilbert space fermion and tensor bases
import numpy as np # generic math functions
##### define model parameters #####
L=4 # system size
J=1.0 # hopping
U=np.sqrt(2.0) # interaction
mu=0.0 # chemical potential

"""
Explanation: The Fermi-Hubbard Model
This notebook shows how to use the tensor_basis constructor to build the Hamiltonian of interacting spinful fermions in 1d, desctibed by the Fermi-Hubbard model (FHM):
$$H = -J\sum_{i=0,\sigma}^{L-1} \left(c^\dagger_{i\sigma}c_{i+1,\sigma} - c_{i\sigma}c^\dagger_{i+1,\sigma}\right) - \mu\sum_{i=0,\sigma}^{L-1} n_{i\sigma}  +U\sum_{i=0}^{L-1} n_{i\uparrow }n_{i\downarrow } $$
where $J$ is the hopping matrix element, $\mu$: the chemical potential, and $U$ -- the onsite $s$-wave interaction. 
We begin by loading the libraries and defining the model parameters:
End of explanation
"""


# define boson basis with 3 states per site L bosons in the lattice
N_up = L//2 + L % 2 # number of fermions with spin up
N_down = L//2 # number of fermions with spin down
basis_up=spinless_fermion_basis_1d(L,Nf=N_up)
basis_down=spinless_fermion_basis_1d(L,Nf=N_down)
basis = tensor_basis(basis_up,basis_down) # spinful fermions
print(basis)

"""
Explanation: To build the basis for spinful fermions, we take two copies of the basis for spinless fermions and tensor them using the tensor_basis constructor. While the tensor_basis can be used to tensor any two bases objects, it does not allow for passing symmetries, other than particle number conservation (we are currently working on developing a separate class which will allow using all symmetries for spinful fermions). 
To this end, we define the number of spin-up and spin-down fermions, and proceed as follows:
End of explanation
"""


from quspin.basis import spinful_fermion_basis_1d

basis = spinful_fermion_basis_1d(L,Nf=(N_up,N_down))
print(basis)

"""
Explanation: Alternatively, one can use the spinful_fermion_basis_1d class as well. This class, unlike the tensor_basis class can handle various 1d symmetries in the usual way and should be preferred for dealing with the FHM.
End of explanation
"""


# define site-coupling lists
hop_right=[[-J,i,(i+1)%L] for i in range(L)] #PBC
hop_left= [[+J,i,(i+1)%L] for i in range(L)] #PBC 
pot=[[-mu,i] for i in range(L)] # -\mu \sum_j n_{j \sigma}
interact=[[U,i,i] for i in range(L)] # U/2 \sum_j n_{j,up} n_{j,down}
# define static and dynamic lists
static=[
        ['+-|',hop_left],  # up hops left
        ['-+|',hop_right], # up hops right
        ['|+-',hop_left],  # down hops left
        ['|-+',hop_right], # down hops right
        ['n|',pot],        # up on-site potention
        ['|n',pot],        # down on-site potention
        ['n|n',interact]   # up-down interaction
                                ]
dynamic=[]
# build Hamiltonian
no_checks = dict(check_pcon=False,check_symm=False,check_herm=False)
H=hamiltonian(static,dynamic,basis=basis,dtype=np.float64,**no_checks)

"""
Explanation: Defining the site-coupling lists is the same as before (mind the signs in the fermion hopping operator, though!). 
The tensor_basis accepts extended operator strings. The idea is that within the subspace of each basis, we use the operator strings belonging to the corresponding underlying basis (for spinless_fermion_basis_1d, the allowed operators are "+", "-", "n", and "I"). We then use a ...|... to separate the operators that act on spin-up (left) and spin-down (right). 
For instance, the hopping operators $c_{j,\uparrow}c^\dagger_{j+1,\uparrow}$ and $c_{j,\downarrow}c^\dagger_{j+1,\downarrow}$ are represented as '-+|I' and 'I|-+', repsectively, where 'I' stands for the identity (and can be dropped, see below). On the other hand, the spin-flip hopping process $c_{j,\uparrow}c^\dagger_{j+1,\downarrow}$ would mix the spin-up and spin-down sectors and would take the form '-|+'.
End of explanation
"""


