I am trying to define a Hamiltonian on a quasi-1D spin-\frac{1}{2} system and am a little confused about the syntax to be used. I found an example for the Heisenberg chain on the documentation site:
lattice = square_lattice(Nx, Ny; yperiodic = true)
os = OpSum()
for b in lattice
os .+= 0.5, “S+”, b.s1, “S-”, b.s2
os .+= 0.5, “S-”, b.s1, “S+”, b.s2
os .+= “Sz”, b.s1, “Sz”, b.s2
end
Is it possible to explain how exactly this notation works? For example, how would one define the following local Hamiltonian on a square lattice
Hi Arnab, I’m not sure from the expression you gave how your Hamiltonian is defined. Are i,j coordinates on a two-dimensional lattice? If so, what is the lattice? Is there a sum over i also, or only over j? Thanks.
I apologize for not being clear. The Hamiltonian is defined on a 2D square lattice and yes the sum is over both i,j which are the coordinates of the sites. I have edited the main post accordingly.
You can think of each element of the lattice object (b in your code above) as a unique bond of the square lattice. In DMRG, we always map the lattice positions to integers since DMRG is ultimately a 1D algorithm, so if you print out b.s1 and b.s2 you will see them running over the integer labels of each site of the square lattice (i.e. something between 1 and Nx * Ny). There is an internal convention for how we map each site of the square lattice to a linear ordering, basically:
Note that you can always print things in your code, for example @show b inside the for loop or type os or lattice at the Julia REPL, which will show you what those objects are and you should see terms with bonds based on the convention I show above.
That helps a lot! The only thing that I am still confused about how to actually define more than two local operator using this notation. I guess this confusion comes from not understanding how to use the .s1 and .s2 attributes properly.
For example, can I ask what the syntax will be for the operator Z_3X_5X_8Z_7 ( in general Z_iX_{i+N_y-1}X_{I+2N_y}Z_{i+2N_y-1} in the 3\times 3 lattice that you mentioned above.
Hi Arnab,
Matt’s drawing above of a 2d system mapped to a 1d system (the numbers define the 1d ordering) could help you here. I can see why for more than two-site operators the b.s1 and b.s2 variables no longer help.
What you would need to do for e.g. the specific lattice Matt wrote is, if you had a operator in 2D notation that was, for example Z_{i,j} X_{i,j+1} Z_{i,j+2} then in the 1D notation it would be Z_{1,1}, X_{1,2} Z_{1,3}. (Here I am making up a 3-site operator that goes up the first column.
So you would just need to work out on paper each of the cases of your operators, and how they map into a “zig zag” wrapping of a 1D system onto a 2D lattice.
You might end up with certain “rewriting rules” similar to the example I wrote above. Like if an operator string goes up a vertical column, the 1D index just goes up by one. But then if it goes horizontally over by one column, then the 1D site index goes up by a value N_y where N_y is the size of the system in the y direction.
So, a rough guess (please check) for your operator at the very first post above could be: Z_{i,j-1} Z_{i,j+1} X_{i,j} Z_{i-1,j} Z_{i+1,j} \rightarrow Z_{i-1} Z_{i+1} X_{i} Z_{i-N_y} Z_{i+N_y}
One other thing to keep in mind is that the OpSum system does not “know” what a 2D system is. All it knows about are 1D systems and when it gets a 2D Hamiltonian it just thinks the user is making a very complicated, longer-ranged 1D Hamiltonian. But in disguise it is a 2D Hamiltonian.
Also note that some of these things will get a bit easier with designs we are working out in GitHub - mtfishman/ITensorNetworks.jl: A package with general tools for working with higher-dimensional tensor networks based on ITensor., where you can make an OpSum with more general site labels (i.e. you will be able to use a 2D grid labelling, and define interactions based on more general graph objects that generalize the lattice objects like square_lattice of ITensors.jl). That package isn’t quite ready for general usage so for now you are better off trying it out how Miles described.
This makes a lot of sense! So, just to confirm, for more than two-site operators, am I better of just ignoring the “lattice = square_lattice(Nx, Ny; yperiodic = true)” construction and just mapping everything back to a standard 1D lattice with long range Hamiltonians?
I see why you would ask that, but the lattice = square_lattice(...) code is in fact also just mapping everything back to a standard 1D lattice with long range Hamiltonians. It is just helping you to do that for the case of nearest-neighbor 2D interactions on a square lattice mapping back to 1D. All Hamiltonians for DMRG are 1D lattice Hamiltonians, possibly long-ranged (if the underlying physics is that of a 2D system).
Hopefully that clears things up for you, but if not, I would encourage you to print out the values b.s1 and b.s2 for the bond (b) objects you get when iterating over the lattice array returned by square_lattice, then compare these values (which are 1De values) to the number in Matt’s drawing above, and you will see they are the same numbers.