Postlifts: mapping scalar mesh vertices to arbitrary hamiltonian parameters

[pablosanjose]

Closes #75

This PR does various changes to meshes and bandstructures

• Disallows non-numeric vertices (i.e. they cannot directly represent non-scalar parameters of ParametricHamiltonians)
• Allows to specify a lift kwarg in bandstructure(h, meshspec; lift = ...) that is now treated as a "post-lift": first, the vertices of mesh = buildmesh(meshspec, h) are lifted to the "embedding mesh", and then the user-provided lift is applied on the result.
• Remove the silent lifting of, say, a 2D marchingmesh to a 3D Brillouin zone. We now require the embedding mesh dimensions to match the hamiltonian's parameter/bloch phase space dimension (EDIT: if no postlift is provided thought the lift option)
• The mesh metric used for finite differences in codiagonalization (for degeneracy resolution) is now performed on the mesh itself (without lifting), as postlifts are not guaranteed to produce scalar parameters with a well defined norm

The "embedding mesh" of a linearmesh((0,0), (0,1)) is a 2D mesh (because vertices are 2D). For a marchingmesh((0,1), (0,1), (0,1)) it is 3D, because there are three axes. So lifting is done in two steps: first from the "cut" to the "embedding mesh" and then however the user specifies with his/her lift. In practice, the user-provided lift (called postlift only in the code, the kwarg is still lift) is the actual mapping from the embedding mesh to the parameter space of a ParametricHamiltonian. To understand this better, here is an example.

Take a parametric graphene Hamiltonian, with a vectorial parameter E that represents an electric field (producing a lattice-periodic potential E' * r)

ph = HamiltonianPresets.graphene() |> parametric(@onsite!((o, r; E) -> o + E' * r))

Now, we can build a bandstructure along a cut in "electric-field-space" like this

bandstructure(ph, linearmesh((0,0), (1,1)), lift = (Ex, Ey) -> (SA[Ex,Ey], 0, 0))

This sweeps E from SA[0,0] to SA[1,1]. First, the linear mesh coordinates are mapped to a 2D mesh space (where the nodes (0,0) and (1,1) live). Then the postlift (Ex, Ey) -> (SA[Ex,Ey], 0, 0) converts the scalar 2D mesh vertices to a static vector, while fixes the bloch phases to zero. The result of the user-provided lift is then always between the embedding mesh and the parameter/bloch phase space.

[pablosanjose]

This required some subtle fixes to get CI back to green.

applylift was forced to produce a SVector before. Since the lift maps onto a parameter+phase space, where parameters are of arbitrary type, that no longer made sense. applylift now produces a tuple, that is passed to bloch! directly. It is then bloch! who does the work of splitting phases from parameters, and upgrading the former to an SVector.

Since here we no longer automatically extend the dimensions of a symbolic linearmesh node such as :Γ to the dimension of the Hamiltonian (since we separate the concept of embedding mesh dimension from the lifted parameter+phase space, with a postlift connecting them), it became impossible to define the linearmesh of purely symbolic nodes independently of h. In the last commit I've changed LinearMeshSpec to allow storing unresolved Symbols for the nodes, and a resolve(::LinearMeshSpec, h) to resolve the symbols back to a mesh vertex of proper dimensions, matching h. The rationale here is that if a symbol is given for a node, we are actually assuming the embedding mesh equals the Brillouin zone.

In any case, if a mix of symbols and tuples are given for linearmesh nodes, they must now be consistent, and symbols are expanded to the dimension of the non-symbolic nodes.

Finite differences in codiagonalizer needed a fix to disallow shifts along parameters when they are non-numeric.

I also added some more useful error message in case a dimension mismatch is found between mesh vertices and Hamiltonian dimensions.

[codecov-commenter]

Codecov Report

Merging #76 into master will increase coverage by 0.10%.
The diff coverage is 90.00%.

@@            Coverage Diff             @@
##           master      #76      +/-   ##
==========================================
+ Coverage   62.31%   62.41%   +0.10%     
==========================================
  Files          15       15              
  Lines        2351     2360       +9     
==========================================
+ Hits         1465     1473       +8     
- Misses        886      887       +1     

Impacted Files	Coverage Δ
src/parametric.jl	82.95% <ø> (+0.19%)	⬆️
src/mesh.jl	91.08% <87.23%> (-2.84%)	⬇️
src/bandstructure.jl	92.30% <100.00%> (ø)
src/diagonalizer.jl	52.38% <100.00%> (-0.38%)	⬇️
src/hamiltonian.jl	68.58% <100.00%> (+0.23%)	⬆️
src/tools.jl	61.97% <100.00%> (-1.05%)	⬇️
src/model.jl	71.72% <0.00%> (+0.52%)	⬆️
... and 2 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 0ef1fb2...01f3dd7. Read the comment docs.

[pablosanjose]

I do realize this PR might be increasing the complexity of our API quite a bit, since the user needs to understand two types of lifts instead of one (i.e. cut -> embedding mesh -> parameter/phase space). The need of this would be even clearer if we had more than one "cut"-like mesh (higher-order generalizations of linearmesh). Hopefully the idea is sufficiently clear, however. I don't see any other way to solve the problem of merging arbitrary non-scalar parameters (a very real usecase) and generalized bandstructures (i.e. keeping the proposal by @BacAmorim of having bloch! operate on a fused parameter/phase space, see #28).

Feedback welcome on the clarity of the change, or otherwise.

[BacAmorim]

So let us assume that a parametric Hamiltonian depends on a matrix. The user wants to plot some band structure along a cut in the bloch/parameter space. So he has to build a mesh with a lift encoding how the mesh points are mapped to the bloch/parameter space (written as a SVector). Then he has do define another lift function that transforms the SVector into a bloch vector and a matrix that parametrizes the Hamiltonian.

If this is true, the user always has to define how a collection of Number parameters must be encoded into a matrix that feeds the Hamiltonian. Why not demand form the start that the arguments of a parametric Hamiltonian must be Numbers? It seems that the complexity of the code increases with no actual benefit. I might be missing something here?

Going to your example. In this proposal, the user does

ph = HamiltonianPresets.graphene() |> parametric(@onsite!((o, r; E) -> o + E' * r))
bandstructure(ph, linearmesh((0,0), (1,1)), lift = (Ex, Ey) -> (SA[Ex,Ey], 0, 0))

Is this better than forcing the use to define the parametric Hamiltonian as:

ph = HamiltonianPresets.graphene() |> parametric(@onsite!((o, r; Ex, Ey) -> o + SA[Ex,Ey]' * r))

?

I understand that in some cases, if might be very cumbersome to write something of the form of (o, r; Ex, Ey) -> o + SA[Ex,Ey]' * r. But in that case the user can define a helper function
electric(Ex, Ey) = SA[Ex,Ey] can then call (o, r; Ex, Ey) -> o + electric(Ex, Ey)' * r. This would essentially be a user implemented postlift. Why do I think this is better:

1. It removes one additional complication to the interface that the user must learn (postlift)
2. Given the constrain "the parameters of a ParametricHamiltonian must be Numbers" and the fact that the argument of modifiers @onsite! and @hopping! are just anonymous functions, the user can easilly come up with a solution suitable for his particular case to make his coding simple to read and write (possibly breaking the function into several ones...).

[pablosanjose]

I did consider this, but I felt that in practice (because I'm actually using this stuff for a project), the verbosity of forcing scalar parameters was suboptimal, since typically the models can be complex. Allowing to write them in terms of arbitrary parameters seemed good.

In any case I think I did not clearly explain myself. The used never has to define two lifts, only one, from the embedding mesh to the parameter/Bloch space. It is only internally that we need to define another lift in cases where the MeshSpec corresponds to a cut of the embedding mesh. But this is transparent for the user. She doesn't need to understand the postlift concept. She only needs to be aware that the source space is the embedding mesh (so it should be clear what that is, given a MeshSpec), and that the target space is parameter/Bloch space, so it is not so complicated. In usual full-bandstructure computations the mapping between these is the identity, so no need to define a lift at all (even when using linearmesh). It is only when using ParametricHamiltonians (or when wanting to implement a funky transformation of the Brillouin zone yourself) that you need to specify a lift.

Does that change your perspective?

[pablosanjose]

One more thought that may clarify things further. Remember you once proposed to merge the lift into the Mesh itself, and I pushed back because I argued that the lift depends on the Hamiltonian? With the current lift/postlift perspective things are clearer. The first lift (from the cut to the embedding mesh) should indeed be part of the mesh object. It is only the postlift that should be separate. So an improvement here would be to refactor the MeshSpec type to be a general "cut mesh + lift to embedding mesh" object. A marchingmesh and linearmesh would both be specific instances of this, both with a clear concept of embedding mesh encoded into the type. Internally one could explicitly store the coordinates of the cut vertices both in the cut manifold and the embedding manifold, or just one of them plus a "prelift" function that maps between them. That should be a transparent implementation detail. Then, we are left with only one (post)lift outside of the mesh object. Maybe that's cleaner.

[pablosanjose]

The latter idea about storing prelifts inside MeshSpecs can be implemented orthogonally to this PR. I will merge this in a bit if there is no further feedback.

[BacAmorim]

Sorry for the late reply (busy week).
I think I still didn't fullly understand the proposal. Let us go to your example of a 2D, with an in-plane electric field. The parameter/bloch space is then (Ex, Ey, kx, ky). Now assume, we want to see the bandstructure along the line: [(x, x, x, x) for x in 0:0.1:1] So we have a 1D coordinate space. And we have a prelift function that goes from the coordinate to the embeeding space:

prelift: Tuple{Float64} -> Tuple{Float64, Float64, Float64, Float64} [math/julia frankenstein notation]
(x) -> (x, x, x, x)

Now, the electric field is actuallly passed as a SVector. So we need a poslift that goes from the embeeding space to the target space of the parametric Hamiltonian:

postlift: Tuple{Float64, Float64, Float64, Float64} -> Tuple{SVector{2, Float64}, Float64, Float64}
(Ex, Ey, kx, ky) -> (SVector(Ex, Ey), kx, ky)

I initially understood that with this proposal the user only actually needs to specify a function that is a composition of both: lift = postlift(prelift), which goes directly from the coordinate to the parametric Hamiltonian target space:

lift: Tuple{Float64} -> Tuple{SVector{2, Float64}, Float64, Float64}
(x) -> (SVector(x, x), x, x)

Is this correct? If so, why would we still want a lift included in the MeshSpec?

[pablosanjose]

Hey @BacAmorim, thanks for your post.

That is not quite correct. The linearmesh((0,0), (1,1)) lives in a 2D embedding space (because the nodes are 2D). At this level there is no connection to the Hamiltonian whatsoever, so it doesn't make sense to say (yet) if this is (kx, ky), (Ex, Ey) or something else. It is a mesh that exists independently of any Hamiltonian. However, it already incorporates a mapping from a line cut coordinate (the distance along the (0,0) -> (1,1) path) and the 2D embedding space. That is the prelift: NTuple{1, Float64} -> NTuple{2, Float64} . This can still be defined independently of the Hamiltonian^{*}.

Now, when we want to "evaluate" the ParametricHamiltonian across this linear mesh, we need another mapping, from the 2D space to the parameter/bloch space. In this case it is postlift: NTuple{2, Float64} -> Tuple{SVector{2,Float64}, Float64, Float64}.

So what this PR does is to first apply the prelift, as obtained by buildlift(meshspec[, h]) (where the h is still needed because of the ^{*} below), and subsequently apply the lift entered by the user, if present (which in this discussion we're calling the postlift, but we don't expose such a confusing name to the user).

---

^{*}: That is true as long as we don't begin to write things like linearmesh(:Γ, :K; samelength = false), in which case we still need some information about the h Hamiltonian Brillouin zone (at least its Bravais vectors). This problem makes the current PR less clean than it should be (we still use h to compute the prelift). Not sure yet how to completely separate prelift and postlift and still keep the symbolic BZ names for convenience.

[BacAmorim]

Ok, I think I have a better understanding of the proposal. But the prelift/postlift destinction (even if not directly exposed to the user) still seems unnecessary to me. For this particular example, it would seem more natural for the user to just specify a linearmesh(0, 1) and then provide a lift = x -> (SVector(x[1], x[1]), x[1], x[1]).

Maybe my confusion is due to the name lift as I would expect this to map a lower dimensional space to a higher dimensional one. With the present proposal, this might not be case, as we can have just a repackaging of the params/bloch phases to a form suitable to the parametric Hamiltonian. For example, going back to the case of 2D system with electric field. We want to obtain the bandstructure in the 4D (Ex, Ey, kx, ky) space. We would use marchingmesh((0, 1), (0, 1), (0, 1), (0, 1)) to create a 4D mesh. With this proposal we would then write lift = (Ex, Ey, kx, ky) -> (SVector(Ex, Ey), kx, ky), so a map between two 4D spaces. So maybe a more descriptive name would be parametrize = (Ex, Ey, kx, ky) -> (SVector(Ex, Ey), kx, ky) or mesh2params = (Ex, Ey, kx, ky) -> (SVector(Ex, Ey), kx, ky) ?

[pablosanjose]

As usual, you make excellent points, so I'm going to delay merging this, and give us some time to elaborate.

In particular, you are implicitly proposing an alternative to the whole cut/lift approach. You say, have a single kind of mesh, parametrized by dimension (we could just make do with marchingmesh, and call it simply mesh(...; type = :marching). We would then implement cuts fully through a single parametrize or mapping or mesh2params kwarg interface (I personally like the name mapping). Then, all bandcut functionality is concentrated on the mapping, and Issue #75 is solved through a single mapping function. So, following your reasoning, these could be some usage examples

Full band, with an L-dimensional h::Hamiltonian

bandstructure(h)
bandstructure(h, mesh((0, 2pi),... ::Vararg{Tuple,L}; points = ...)) 

Straight bandcut, with an L-dimensional h::Hamiltonian

bandstructure(h, mesh(0,1), mapping = k -> (k, k,... )::NTuple{L})

Piecewise bandcut along a 1D cut from k = (0,0) to k = (0,pi), with an 2-dimensional h::Hamiltonian

bandstructure(h, mesh(0,1), mapping = piecewise((0,0), (0,pi)))

Piecewise bandcut along a 1D cut with symbolic nodes

bandstructure(h, mesh(0,2), mapping = piecewise(:Γ, :K, (pi,pi)))

General parametric Hamiltonian ph, hand-written mapping

bandstructure(ph, mesh((0,1), (0,2)), mapping = (x, y) -> (SA[x,x+y], y, y^2))

Parametric mapping from a 1D mesh, with the ph(; E) in the OP, and a sweep from E = SA[0,0] to E = SA[1,1]

bandstructure(ph, mesh(0,2), mapping = piecewise((SA[0,0], 0, 0), (SA[1, 1], 0, 0)))

The piecewise function could produce a function that linearly interpolates between arguments (using broadcast .+ and .*), as the arguments. If given purely symbolic nodes, it could build a lazy object that only resolves into a concrete function when the h becomes known.

Please let me know how you feel about this alternative. It does indeed seem simpler... If you agree, I might close this and start a new PR where we can polish the details.

[BacAmorim]

I really like this! I also like the name mapping!
I think that this is single unifying API that appears to be quite natural and solves several problems in one go. I imagine that someone that is trying to learn Quantica (besides me) will have no troubles understading this feature.

One question: with the piecewise mapping, now are you going to split the points of the 1D mesh through the different segments? As equally as possible between the different segments, based on the length of each segments? (While I like this and have proposed it before, this would require defining a metric for arbitrary parameters, which was one of the initial problems here). linearmesh solved this by allowing explicitly specify the number of points per segment. How would this be addressed here? Maybe keep linearmesh and marchingmesh? Also, with the piecewise option, there is really no difference in writing

mesh(0,1), mapping = piecewise((SA[0,0], 0, 0), (SA[1, 1], 0, 0))

vs

mesh(0,2), mapping = piecewise((SA[0,0], 0, 0), (SA[1, 1], 0, 0))

The only thing that matters is the number of points in the segment. Correct?

I would suggest keeping in the API both linearmesh and marchingmesh (maybe just calling it mesh as the the marching part is an implementation detail). But I I am sure you will come up with a more elegant solution.

[pablosanjose]

Ok, let's move this over to an independent issue