## Lattices and Hamiltonians

NOTE: this is very out of date!

The new structure (from version 0.7.4 onwards) uses a recursive tree structure to define semi-periodic lattices and operators.
The basic operations are `repeat`

, to repeat some lattice `N`

times, and `join`

, to join two dissimilar lattices.

For example, the {$SO(4)$} symmetric Hubbard model is defined on a bipartite lattice, so we need to join two single-sites together into a unit cell.

#include "matrixproduct/lattice.h"

#include "matrixproduct/mpoperatorlist.h"

#include "matrixproduct/operatoratsite.h"

#include "models/hubbard-so4.h"

int main()

{

Lattice SiteA = CreateSO4HubbardSiteA();

Lattice SiteB = CreateSO4HubbardSiteB();

Lattice UnitCell = join(SiteA, SiteB);

The `CreateSO4HubbardSiteA()`

and `CreateSO4HubbardSiteB()`

functions return objects of type `SiteBlock`

. There is a conversion constructor to construct a Lattice object from a single `SiteBlock`

. (As a convenience, there are constructors to make the `join`

of multiple `SiteBlock`

's too - see `models/hubbard-so4.cpp`

for an example.) If we want to supply a `SymmetryList`

object, we can do that here too, at the point of constructing a `Lattice`

out of one or more `SiteBlock`

's.

We now want to repeat the unit cell {$N=L/2$} times, to construct the full lattice:

int L = 50; // the lattice size

Lattice MyLattice = repeat(UnitCell, L/2);

if (L%2 == 1)

MyLattice = join(MyLattice, Lattice(CreateSO4HubbardSiteA()));

Note the trick at the end: if the lattice size is odd, we join a single site (half a unit cell) onto the end of the lattice.

The `repeat`

and `join`

functions can be nested arbitrarily.

Given this lattice file, we can construct an `OperatorList`

,

OperatorList OpList(MyLattice);

OperatorAtSite<OperatorList const, int> C(OpList, "C");

OperatorAtSite<OperatorList const, int> CH(OpList, "CH");

OperatorAtSite<OperatorList const, int> P(OpList, "P");

MPOperator& Hamiltonian = OpList["H"];

There is some confusing naming here; the `OperatorList`

is what is actually stored on the lattice file on the disk. This does contain the `Lattice`

though, available through the `OperatorList::GetLattice()`

function.

The `OperatorAtSite`

objects are convenience functions. With these definitions, The expression `C(5)`

is equivalent to `OpList["C(5)"]`

. The big advantage is that, if we want to replace the '5' with a variable, we can use it directly in an `OperatorAtSite`

, as in `C(x)`

, whereas if we use the `OperatorList`

directly, we need to convert the site label into a string.

With the new lattice system, we barely need the `OperatorList`

at all. In fact, rather than looping over all sites, we can construct the Hamiltonian as

MPOperator Hopping = -t * CreateRepeatedOperator(MyLattice, "CP", "CH");

MPOperator Coulomb = 0.25 * U * CreateRepeatedOperator(MyLattice, "P");

Hamiltonian = Hopping + Coulomb;

The `CreateRepeatedOperator()`

functions are new in version 0.7.4; these take advantage of the lattice structure to construct an optimal periodic (open boundary condition) representation of the operator, repeated on every site.

The efficiency of `CreateRepeatedOperator()`

depends on having a periodic structure in the lattice. If, instead of using the `repeat()`

function, we `join`

'ed the unit cell together {$N$} times, the `CreateRepeatedOperator`

function would need to loop over the {$N$} unit cells, and also use {$O(N)$} memory.

Unfortunately, `CreateRepeatedOperator()`

does not yet utilize the Fermionic commutation relations built into the site operators, so we need to explicitly include any needed {$P \equiv (-1)^N$} terms as required. Hence, we have specified the hopping term as `("CP", "CH")`

, using the operator `CP = C . P`

, which includes the fermion anticommutation operator on the first site. Hence, it is now important to be aware of the ordering of the lattice sites. For non-Abelian symmetries, the `CreateRepeatedOperator()`

function includes the necessary coefficient to construct a dot product, so we do not need the factor 2 that was used previously.

Constructing operators by looping over sites still works, however. We can use this to construct bond operators for the Suzuki-Trotter decomposition. For example,

OperatorAtSite<OperatorList, int> Bond(OpList, "Bond");

for (int i = 1; i < L; ++i)

{

Bond(i) = -t * dot(C(i), CH(i%L+1)) + (U*0.125)*P(i) + (U*0.125)*P(i+1);

}

Bond(1) += (U*0.125)*P(1);

Bond(L-1) += (U*0.125)*P(L);

Here, we have used the new `dot()`

operator, which is the dot product of the two operators, and we have evenly distributed the Coulomb terms over the bond operators. For abelian symmetries, this is equivalent to the ordinary product. For non-abelian operators, this takes care of the coefficients that differ from taking the scalar coupling versus the dot product.

### Future

It is intended that more functions similar to `CreateRepeatedOperator()`

will be added. For example, it would be useful to be able to construct the sum of all even (or odd) bond terms directly, similarly to how we constructed the Hamiltonian. Also, the new lattice format no longer has a provision for naming the lattice sites. This will be fixed very soon - the only question to resolve is what syntax should be used.

The syntax for constructing lattices containing Y-junctions also needs finalizing. This needs to be something like `junction(LeftLead, TopLead, RightLead)`

. Or maybe something like `junction(LeftLead).with(TopLead, RightLead)`

? There is an asymmetry, in that the actual junction will be the right-hand edge of one chain, with the left-hand edges of two other chains. The notation for constructing the lattice should reflect this.