Introducing the current Operator

[pablosanjose]

There is a long-standing design issue around how to represent quantum operators that are not Hamiltonians. We could opt to rename Hamiltonian/ParametricHamiltonian/AbstractHamiltonian to Operator/ParametricOperator/AbstractOperator, and also rename hamiltonian to operator. However that seems a design mistake. All the language we use for modelling (onsite, hoppings...) are actually meant to represent Hamiltonians, and that is actually what we want to model most of the time. It would also be breaking.

We could also just rename the structs, and keep the hamiltonian API unchanged. The output would just be an AbstractOperator, without any type information indicating it is actually a Hamiltonian. But then it would be weird to have GreenFunctions accept any generic operator that is not an AbstractHamiltonian.

In this PR we take a more conservative approach: we keep the AbstractHamiltonian as the central object, and make Operator a derived object that is constructed from an AbstractHamiltonian.

struct Operator{H<:AbstractHamiltonian}
   h::H
end

This is completely non-breaking and keeps all the present structures and APIs untouched, so we can experiment for a while with this, and revert if we find the approach lacking.

As a first example we introduce an API method current(h::AbstractHamiltonian; charge = -I, direction = 1), which returns an Operator representing the current Jᵢ = ∂h/∂Aᵢ along a given direction = i.

CC @BacAmorim

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests ✅

Comparison is base (1d9f7f1) 92.60% compared to head (3336ce9) 92.36%.

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #227      +/-   ##
==========================================
- Coverage   92.60%   92.36%   -0.24%     
==========================================
  Files          34       34              
  Lines        5460     5477      +17     
==========================================
+ Hits         5056     5059       +3     
- Misses        404      418      +14     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[BacAmorim]

Keeping the Hamiltonian as the central quantity and defining operators on "top of" the hamiltonian makes sense to me! This replicates the way we typical model system: define a Hamiltonian and they compute observables from it.

I assume that this would also work for ParametricHamiltonian, correct?

Could you please explain the charge = - I keyword argument?

[pablosanjose]

I assume that this would also work for ParametricHamiltonian, correct?

Yes, absolutely, since ParametricHamiltonian <: AbstractHamiltonian. So you can do J = current(h::ParametricHamiltonian), and then J(blochphases; params...) and you get your sparse current matrix. You can also do Quantica.call!(J, blochphases; params...) to do the same thing in-place, without allocating a new matrix.

Could you please explain the charge = - I keyword argument?

Ah, that's just to allow you to compute generalized currents, such as spin currents. I'm not sure I did the generalization correctly though, I will need to check again.

[pablosanjose]

When writing tests for this PR, I realised that getindex(::ParametricHamiltonian,...) invariably calls getindex on its (non-parametric) parent, which leads to surprising results. This is particularly true with Operator, such as a current. I changed the method to getindex(::ParametricHamiltonian,...; params...), which now applies params, and then calls getindex on the result.

And with that, this is good to go.

[BacAmorim]

The construction of the current operator is simply a modification of the original Hamiltonian hoppings and onsites terms. The structure of more general operators is completelly unrelated to the Hamiltonian (for example a dipole operator).

For this case, it would be useuful to be able to to this:

h = (some Hamiltonian)
operator_model = onsite(...) + hopping(...)
op = operator(h, operator_model) #this operator inherits the lattice from the hamiltonian h, but with onsite and hopping terms defined from the operator_model 

Is this possible at the present?