Link-free redesign - Githubissues

This is a comprehensive redesign of the Elsa internals and API. It addresses some of the comments in #3.

The problem(s)

The design in current master is based in a clear separation between linking of sites (as a geometric property) and the hoppings that are applied to said links (as a physical property). This has a number of drawbacks that became clear by considering the comments in #3 by @Datseris.

One of them is the difficulty to specify different hoppings at different links that differ on, say, a given intercell distance. The only way was to make a function of two positions that would include a condition on said distance, but it's not the most elegant approach.
Another problem is the redundancy between the stored links (of type Links, innternal to Lattice) and the Hamiltonian in System. These contain quite redundant information, although for some problems (say, deforming a lattice), the information in Links was easier to manipulate.
Apart from the redundancy, and the associated memory footprint, there is the performance of building a large system. Both the Links and the Hamiltonian had to be constructed, which leaves at least a factor of ~2 in performance on the table.
One more issue raised, which is debatable but I think I agree with, is that the use of types for all possible directives is not the most idiomatic/Julian approach. It had the advantage of allowing building a system in a way that depended on the applied directives. However, there is a much more Julian way, which is to chain a series of operations on a given system/lattice. This is the approach that this PR now explores

The proposed solution

This PR excises all the Links machinery, and stores all the information of the lattice connectivity directly as hoppings in a Hamiltonian. As such, linking two sites is just a matter of adding a hopping between them to the Hamiltonian. This is similar to the Kwant and Pybinding approach, I believe.

As a result there is no longer a distinction between Lattice and System. So Lattice is gone. A System is just a collection of sublattices, a Bravais matrix and a Hamiltonian (of type Operator). Systems are build with simple constructors that take either 1- Sublats, Bravais vectors and a Model 2- A preset, and optionally a Model to be applied to the preset 3- Another System and a Model to replace the Hamiltonian of the original system.

In addition to this constructor, the System API includes a number of functions, that can be applied sequentially in any given order to build derived systems from one or more systems. They allow piping with |> for cleaner syntax. These functions are modifysublats!, grow, transform, transform! and combine (for the moment). We also have the function hamiltonian to extract the Hamiltonian of a system. They are all documented, both with docstrings and in the new Getting Started doc.

One problem I had with this PR is that it required almost a full rewrite of the Elsa basics, and as a result the work is not yet at feature-parity with master. This PR allows to build Hamiltonians, and to visualise systems with the new PR in ElsaPlots.jl, but bandstructure calculations (to be renamed "Spectrum") is not yet done.

Performance-wise this PR is a step forward. Building a Hamiltonian from scratch is now approximately 1.5x to 2x as fast as in master, and uses less than half as much memory. Regarding the new syntax, I think I like it better than before, although it is less compact, but I would appreciate some feedback if anybody can provide some.

Some runtime/syntax comparisons (circular graphene flake of radius 100a0, nearest-neighbor hopping)

Current master:

julia> @time using Elsa
  0.853706 seconds (1.80 M allocations: 114.473 MiB, 8.07% gc time)

julia> test() = System(Lattice(:honeycomb, LinkRule(1/√3), Region(:circle, 100)), Model(Hopping(1)))

julia> using BenchmarkTools; @btime test()
  36.359 ms (2176 allocations: 43.33 MiB)
System{Float64,2,0} : 0D system in 2D space with Float64 sites.
    Bravais vectors : ()
    Number of sites : 72562
    Sublattice names : (:A, :B)
    Unique Links : 108445
    Model with sublattice site dimensions () (default 1)
    Bloch operator of dimensions (72562, 72562) with 216890 elements

This PR

julia> @time using Elsa
  0.803117 seconds (1.74 M allocations: 111.029 MiB, 9.03% gc time)

julia> test() = System("honeycomb", Model(hopping(1, range = 1/√3))) |> grow(region = Region("circle", 100));

julia> using BenchmarkTools; @btime test()
  20.890 ms (946 allocations: 13.41 MiB)
System{Complex{Float64},Float64,2,0} : 0D system in 2D space
  Bravais vectors     : ()
  Sublattice names    : ("A", "B")
  Sublattice orbitals : (1, 1)
  Total sites         : 72562 [Float64]
  Total hoppings      : 216890 [Complex{Float64}]
  Coordination        : 2.989030070835975

CC @Datseris

Hi Pablo!

This PR excises all the Links machinery, and stores all the information of the lattice connectivity directly as hoppings in a Hamiltonian. As such, linking two sites is just a matter of adding a hopping between them to the Hamiltonian. This is similar to the Kwant and Pybinding approach, I believe.

Yes that is correct.

I read through this PR message. I very happy that you got almost ~2 performance benefit (which for tight binding applicaitons is huge as they "human measurable" time, i.e. more tha miliseconds). Also happy to see that the core types necessary to get going are just 3: namely the sublattices, the Bravais vectors and the Model. This sounds very intuitive to me, and if we strip the process to the bare basics this is the minimum amount of information one needs! it is also nice to see taking advantage of the pipe syntax. Most users coming from Python etc. are unfamiliar with it of course, but I believe it is worth it to take "risks" let's say that may lead to a nicer interface. I personally like your approach even though I never use the pipe syntax personally.

I haven't read the documentation yet. I will do that and will give you more feedback!

Ah wait I am a bit confused. Why is it that the example has:

Bravais vectors     : ()

? Shouldn't it list the actual vectors?

You mean example 2, yes? That is a graphene flake of a given radius, so it is not an infinite Bravais lattice. It has no translational invariance, and no Bravais vectors. It is very different from System("honeycomb"), which is truly an infinite system with two Bravais vectors.

Alright I am going through the docs now! You should consider to answer all the questions I ask here in the documentation document, to not do double effort!

hopping functions take the position r of the center of the link and its vector dr (h = (r,dr) -> hfunc(r,dr))

What does "center of link" mean as r? Is it the position of Site A while dr is the difference between A and B? If yes then the way it is phrased in the docs is unclear. The way I read it in the docs is that r (as a vector) is r = (r_A + r_B)/2 which is certainly unintuitive.

Is it mandatory to use a specific type for the sublattice naming scheme, e.g. only Strings or only Symbols?

Or one can use anything?

The default values of Model should be exposed more, and more clearly I believe. (See below where I ask the value of onsite entries).

Does transform! keep the hoppings as is?

On the first example:

Create another sublattice :B, again with a single site, but at (0.0, 0.5/sqrt(3.0)

Typo here, because you use "B" instead of :B in the code snippet.

model = Model(hopping(1, range = a0/√3));

What is the value of the onsite elements? By default zero? How does the "Model" know when to stop adding onsite elements? This instruction only has info about the hopping. This has confused me a bit. Does Model always correspond to an infinite system that must be made finite afterwards?

From your comments on this PR and from the documentation I think it is clear that now there is no such thing as "lattice". A system is instantiated when and only when the Hamiltonian of the system gets an entry at some position.

In addition the docs say:

Model: A tight-binding model given as a set of hopping and onsite instructions Model(onsite(...), hopping(...), ...), with the following syntax

So I guess this expression:

System("honeycomb") |> grow(region = Region("circle", 3))

Just uses the already existing functions in System to "grow" it? Where does it get the onsite values from?

But still, the most important question is:

As it seems, the first instance of System a user can create is always infinite and you must "grow" it to get a finite system . Correct? If yes, then adding this exact statement in the docs could be useful.

For later analysis, is it possible to still obtain the "original" Bravais vectors from a finite system?

The "presets" part can be improved I think. It is not very good practice to have different presets (all of type string) give different type of output (of course this is totally personal opinion). You mention that some presets do and some do not have a "model". After reading the presets parts, it is unclear to me what a preset is.

Conclusion

This looks seriously good. I actually like the piping syntax. But it should be made more clear what is finite and what is infinite I believe. The name grow is truly opposing its name. grow takes an infinite system and makes it finite. So certainly the system doesn't "grow" in reality :P

And in general I think this:

grow(system; supercell, region) to expand the unit cell of a system with Bravais vectors bravais so that the result has bravais vectors bravais * supercell (could be less vectors than the original). region is a boolean function of spatial position that specifies if a site at r in the unit cell should be included in the result.

Can be a bit confusing for a newcomer. Shouldn't it be simple to state how to make a "finite" system? This involves the concept of Super cells, and the product of a bravais vector with a supercell. I admit that I myself do not immediatelly understand what it means. But I do immediatelly understand "hey this is a finite graphene system looking like a circle!".

lastly, as a personal comment:

julia> grow(sys, region = incircle)
System{Complex{Float64},Float64,2,0} : 0D system in 2D space
  Bravais vectors     : ()

I will never find intuitive the choice of not displaying the Bravais vectors. Yes, on paper you are absolutely true, there are no "true Bravais vectors", but still there are the vectors that connect the unit cells with each other and these are invariant! You already state clearly that the system is 0D, in the first line. Is there no way to persuade you to just normally display the Bravais vectors as well?

I'll be pushing a commit which includes some changes to the docs that address most of these. However here are some specific answers.

What does "center of link" mean as r? Is it the position of Site A while dr is the difference between A and B? If yes then the way it is phrased in the docs is unclear. The way I read it in the docs is that r (as a vector) is r = (r_A + r_B)/2 which is certainly unintuitive.

You read it correctly. It's now clarified in the doc. The choice is very intentional. After extensive use of my previous package MathQ I've found that what you nearly always need is link center and link vector. In particular, a Peierls phase is exp(i dot(A(r), dr)), which gets messy with r1 and r2

Is it mandatory to use a specific type for the sublattice naming scheme, e.g. only Strings or only Symbols? Or one can use anything?

It one or the other. I cannot seem to settle on one. What do you think? I've pushed a commit that changes it back to Symbols, which I think is visually lighter and faster to type. It's probably just a matter of choosing one and sticking to it.

Does transform! keep the hoppings as is?

Yes. For a transform that changes hoppings my plan is to implement a new transformation that's probably going to be called field, and which allow to implement things like magnetic fields etc very fast, without recomputing the Hamiltonian Operator. That's not done yet, however.

What is the value of the onsite elements? By default zero?

Yes, anything not specified in Model is zero. In particular Model() creates a Hamiltonian full of zeroes.

How does the "Model" know when to stop adding onsite elements? This instruction only has info about the hopping. This has confused me a bit.

No, this confusion is probably due to your prior experience with Kwant. In Elsa.jl there is a lattice, regardless of its Model: it is defined as the site positions in the unit cell plus the Bravais vectors. What is probably confusing you here is that when you compute the Hamiltonian of such a lattice you don't create a "big" Hamiltonian corresponding to many unit cells of said lattice. It is a Bloch Hamiltonian, whose dimension is the same as the number of orbitals in a unit cell. It has furthermore a dependence on wavevector in the infinite lattice. Hence the Model fills onsite elements only for the orbitals in the unit cell. The hopping is more tricky, as it applies to the unit cell, and also to links between neighboring unit cells (the latter are precisely what make the Bloch Hamiltonian depend on wavevector).

Does Model always correspond to an infinite system that must be made finite afterwards?

No, see above. The model can be applied to a finite or infinite system, and produces always a generalization of a finite-system Hamiltonian, which is the Bloch Hamiltonian, and which of course coincides with the usual Hamiltonian in finite lattices.

From your comments on this PR and from the documentation I think it is clear that now there is no such thing as "lattice". A system is instantiated when and only when the Hamiltonian of the system gets an entry at some position.

No, this is not correct, but the meaning of lattice is probably not what you have in mind. The lattice, as defined above, exists even when you use a model with all-zero elements, i.e. an empty Model(). The lattice contains the information about the positions and periodicity in space of the lattice (periodicity of both the geometrical positions but also of its Hamiltonian). The Hamiltonian contains the information about the linking and onsite energy. So the approach here is not like in Kwant, where as I understand a site is "created" by adding an onsite energy to it.

If you want to create an arbitrary onsite energy on a number of sites in different unit cells you have two options: enlarge the unit cell using grow(supercell =...) or in the future use a field.

System("honeycomb") |> grow(region = Region("circle", 3))

Just uses the already existing functions in System to "grow" it? Where does it get the onsite values from?

The System("honeycomb") has an empty Model(), so all its hoppings and onsite energies are zero (to override just do System("honeycomb", Model(...)). The grow function just repeats the unit cell in space to fill a specific region, which results in a finite system. If the original System had a non-zero Hamiltonian it will get replicated throughout all the region. Alternatively, one can grow with an empty model (which generates no onsite or hoppings in the finite system) and then apply a Model. There are quite a lot of possibilities using just this simple API.

As it seems, the first instance of System a user can create is always infinite and you must "grow" it to get a finite system . Correct? If yes, then adding this exact statement in the docs could be useful.

No, this is not so. The system is infinite if you give it a non-zero number of Bravais vectors. Hence, if you do System(Sublat((1,2))), and specify no Bravais vectors or Model, your system is just a single site at position (1,2), not repeated in space and with zero onsite energy. If you do System(Sublat((0,0)), Bravais((1,0))) instead you are creating an infinite 1D chain in 2D space along the x direction.

It would be good to understand how the docs could be improved to avoid these misundertandings. Do you have any specific suggestion?

For later analysis, is it possible to still obtain the "original" Bravais vectors from a finite system?

You just keep the original system. You do sys1 = System("honeycomb") and then sys2 = grow(sys1, region =...). sys2 does not contain information about the Bravais vectors of sys1, because in general this will not have any meaning for sys2. I know you are thinking: but hey, sys2 is still locally a honeycomb, right? That is right, but it ceases to be so as soon as you add more manipulations to sys2, say do a transform or add vacancies. So it doesn't make sense to me to keep the original Bravais if the new system can be modified.

The kind of thing that you perhaps actually have in mind, is closer to a future feature that will add a type of system, that is likely going to be called LazySystem (please suggests a better name if you can, I'm not sold on this). This is a non-materialised grow command, namely a periodic system together with a supercell and a region, that can however be used to perform Hamiltonian * vector multiplications without constructing the full Hamiltonian of the grown system (hence LazySystem). The purpose of this object is connected to some future large-scale Kernel polynomial method functionality that is still in design stages.

I will never find intuitive the choice of not displaying the Bravais vectors. Yes, on paper you are absolutely true, there are no "true Bravais vectors", but still there are the vectors that connect the unit cells with each other and these are invariant!

They are not invariant once you grow the system. They will be invariant in a LazySystem (or whatever we call it), so indeed, the latter will remember its original Bravais vectors.

The "presets" part can be improved I think. It is not very good practice to have different presets (all of type string) give different type of output (of course this is totally personal opinion). You mention that some presets do and some do not have a "model". After reading the presets parts, it is unclear to me what a preset is.

Ok, how would you design a library of pre-built systems for convenient use? That's all I tried to do with the presets. So a preset is a prebuilt system, that can optionally take some keyword parameters. The fact that most of them have no model is just because there are just a few, and most of them are just the basic Bravais lattices.

grow takes an infinite system and makes it finite. So certainly the system doesn't "grow" in reality :P

That is just one possible use of grow. This is a deceptively powerful function. It just fills a region, making the system finite, if you don't specify the a superlattice kwarg. The essence of grow is to expand the unit cell, while at the same time defining new Bravais vectors, possibly less than the original lattice or none. Hence the name grow, which I quite like (short and descriptive). Check out the docstring of grow:

grow(system::System{Tv,T,E,L}; supercell = SMatrix{L,0,Int}(), region = r -> true)

  Transform system into another system with a different supercell, so that the new Bravais matrix is br2 = br * supercell, and only sites with region(r)
  == true in the unit cell are included. 

  supercell can be given as an integer matrix, a single integer s (supercell = s * I), a single NTuple{L,Int} (supercell diagonal), or a tuple of
  NTuple{L,Int}s (supercell columns). Note that if the new system dimension L2 is smaller than the original, a bounded region function should be
  provided to determine the extension of the remaining dimensions.

Thanks your answers answer all my questions very nicely. Please give me more time to think about this further especially the grow part. Everything else I believe I have understood pretty clearly now.

No, thank you George, your comments are much appreciated. If you find a way to challenge the design further don't hesitate to post.

When designing the fields machinery I realised that this design (and master as well) imposes some limits to future performance. The reason is that System types do not store static information on the orbital dimensions of each sublattice (provided using the norbitals kwarg in Sublat This info is merely stored in a vector inside system.sublatsdata. I just pushed another large commit to address this problem.

The updated PR doesn't take any norbitals kwarg in Sublat anymore. The number of orbitals is derived from the Model itself, which makes much more sense. Additionally, when building the System with the model, care is taken to store sufficient static information in the System parameters themselves to be make orbital dimensions statically inferable.

The new commit also restores the Lattice type (as a bundle of Sublats and a Bravais matrix), so that System is built with a Lattice and a Model. This has no impact on the API, as building an intermediate Lattice is not necessary for the user (Sublats and Bravais can be passed to the System constructor directly, which takes care of building a Lattice behind the scenes).

It also adds further refinements to the docs and docstrings, and fixes some small bugs. Performance stays roughly the same.

Closes #3

pablosanjose / Elsa.jl

Link-free redesign #6

The problem(s)

The proposed solution

Conclusion