It is actually possible, at least in the Z direction, for a device. See e.g.
http://quantumwise.com/forum/index.php?topic=2172.0Now, the real difficulty is that in a device, the translational invariance is broken, or put in another way we don't have periodic boundary conditions. Therefore, we cannot define a band structure in the Z-direction, and therefore we also cannot define the bottom of the CB or top of the VB. The only information we have is the density of states, which we can use at least to get an idea of where the bands are. However, note that we actually get more information than in the classical picture where the band edges are simply lines separating occupied and unoccupied states. We are in really in a position to compute the low occupation in the band gap, due to thermal broadening and tunneling (the exponential tail of the wave function).
In this way you can get a picture like the one attached, for an InAs FinFET-like device under finite source-drain and gate bias (cf.
http://dx.doi.org/10.1109/SISPAD.2013.6650654). This picture uses the LDOS and the results are based on DFT; you can alternatively try the attached scripts to simulate a n-p-n Si junction with a simple sp3s* tight-binding model. In this case the DeviceDOS is used (a bit faster to run).
I put all in one script, to make it easier to run for you, but this it rather long and a bit odd (reading and writing data to file). Normally I split up the computation, averaging, and plotting as 3 separate scripts. The script produces the second attached plot and takes about 10-15 minutes to run on my laptop. Most of this time is however spent computing the density of states, and would parallelize beautifully on a cluster up to 30-50 cores.