Implemented a bulk-shift of the eigenspectrum based on occupations

We implement a generalized method that allows shifting the band-structure based on the velocity of the states.

The first thing we do is to calculate the velocities of each state

<psi_i | dH - e dS | psi_i>

which needs to be done along each Cartesian direction. Since there may be degeneracies it is vital that the we decouple degenerate states. This is done by calculating the matrix

<psi_j | dH - e dS | psi_i>

and diagonalizing to generate the new psi'_i.

After having the velocities we project the velocities along a direction to find the states projection along a given potential direction. If it is positive we can shift the eigenvalue down by a pre-defined "bias" and thus allow this state to enter the valence band. Additionally states with a negative projection will be shifted up in energy (since they are emptied). This is done on all eigenvalues and finally the density matrix is regenerated based on these re-aligned eigenvalues.

All this is controlled via:

BulkBias.V BulkBias.Direction BulkBias.Tolerance

where the last option defines the tolerance for determining a projection to be positive or negative.

It currently only works for Diag.ParallelOverK (since the systems are typically rather small).

The implementation is currently described in this pre-print (submitted): https://arxiv.org/abs/2202.02358

It is very simple and perhaps some of its functionality might be valuable as general output.

For instance one can calculate the velocities of eigenstates and write these to a file for post-processing. The code is there, and easy to use.