# Tight binding packages for Fermi velocity of slabs

I'm studying some topological materials through first principles methods (very trendy I know) and have a question about extracting the Fermi velocity. I'd like to use a wannierised model to analytically calculate the Fermi velocity of a surface states of a slab, but I can't find any existing packages with this functionality. PythTB has nice slab functionalities for TIs and W90 can tell you the Fermi velocity, but only in the bulk.

It's clear to me that the velocity operator is constructed as

$$\hbar v^\alpha_{n} = [U^\dagger \frac{\partial H}{\partial k_\alpha} U]_{nn}$$

with $$U$$ transforming between Wannier/Bloch basis and the derivative analytically expressible as (https://journals.aps.org/prb/pdf/10.1103/PhysRevB.75.195121)

$$\left[ \frac{\partial H}{\partial k_\alpha}\right] _ {ij} = \sum_{R} \exp(ikR) (iR_\alpha)\langle0i|H|Rj\rangle$$

but I'd like to avoid reinventing the wheel if possible. So if anyone knows any packages that implement both slab models and the velocity operator I'd greatly appreciate it. Does this exist? Is this utterly trivial and I'm just being silly/lazy?

Many thanks

I had more of a look and yeah, none of the available packages do this to my knowledge so I just wrote my own code. It isn't particularly hard, more just tedious efficiently figuring out where in your slab hamiltonian each coupling has to go; as you might expect $$H$$ decomposes nicely into a block structure encoding transitions between layers of the slab. For example, for a slab along z (ie 1x1xn) and letting $$(k_x, k_y) = k_{||}$$, the Hamiltonian looks something like
$$\begin{pmatrix} H_{\Delta z = 0}(k_{||}) & H_{\Delta z = -1}(k_{||}) & H_{\Delta z = -2}(k_{||}) & \ldots \\ H_{\Delta z = +1}(k_{||}) & H_{\Delta z = 0}(k_{||}) & H_{\Delta z = -1}(k_{||}) & \ldots \\ \ldots \end{pmatrix}$$
Edit for a point of clarity: This method does indeed work but make sure to read the paper carefully. The diagonal elements of $$[U^\dagger \frac{\partial H}{\partial k_{a}}U]$$ coincide with the velocities but this object is not the velocity operator, as explained a little opaquely in sections B2 and C1 of the linked paper. The object you calculate will have off diagonal elements, but these are expressible as an entirely off diagonal operator related to the Berry phase.