%matplotlib inline
import os
import macrodensity as md
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

extrema = md.vasp_tools.get_band_extrema('OUTCAR_ZnO')
print(extrema)
extrema = md.vasp_tools.get_band_extrema('OUTCAR_ZnS')
print(extrema)

input_file = 'LOCPOT'
lattice_vector = 3.6
output_file = 'planar.dat'

"""
Explanation: Heterojunction offsets
Heterojunctions are the basis of much of the semiconductor technology. From computer memory, to flat-panel display to photovoltaics, semiconductor heterojuctions underly much of the technology created in the past 100 years. Previosuly we constructed band-alignment-diagrams using the IP and EA from slab calculations. However, the details of the interface can also affect the line-up, from charge transfer and lattice strain effects, for example. For an excellent overview of the subject, I recommend: Band engineering at interfaces: theory and numerical experiments.
Our system
For this demonstration we will choose a well-matched interface system. The interface between wurtzite structured ZnO and ZnS. In this case the lattice constants have been set to be those of ZnO. Obviously, due to deformation potentials the value of the offset depends on the lattice parameter used.
Note The values in this tutorial are from very approximate calculations and no interface relaxation was allowed, therefore the absolute actual numbers should not be treated as meaningful.
<img src="HJ.png">
Procedure
We present two procedures for obtaining the offset from a heterojunction calulcation. 
(i) We will intitally calculate the offset for the interface as we have modelled it, this essentially involves the calculation of the macroscopic average of the potential.
(ii) We will then look at how the effects of strain at the interface can be accounted for to yield a "natural band offset". As outlined in Appl. Phys. Lett. 94, 212109 (2009).
Procedure (i)
$VBO = \epsilon_{vbm}^a - \epsilon_{vbm}^b + \Delta V$
$VBO$ is the offset. $\epsilon_{vbm}$ are the eigenvalues of the highest occupied bands from bulk calculations of the two materials (ZnO and ZnS) and $\Delta V$ is the offset in the potential across the interface.
We obtain the $\epsilon_{vbm}$ values from the bulk OUTCAR files as in the slab model exercise. 
Interface calculation
Now we do a calculation of the interface to get the potential profile. Important settings for the INCAR file:
LVHAR = .TRUE.  # This generates a LOCPOT file with the potential

In your example directory there should already be a LOCPOT file. This is the one we will use to analyse the potential and extract the value of $\Delta V$.
In the sample PlanarAverage.py file, all we have to edit are the top three lines. Of these the only one that is not obvious is the lattice_vector parameter. This is just the periodicity of the slab in the direction normal to the surface. In the picture below, this is just the distance between the layers of ZnO.
<img src="HJ.png">
End of explanation
"""


if os.path.isfile('LOCPOT'):
    print('LOCPOT already exists')
else:
    os.system('bunzip2 LOCPOT.bz2')

vasp_pot, NGX, NGY, NGZ, Lattice = md.read_vasp_density(input_file)
vector_a,vector_b,vector_c,av,bv,cv = md.matrix_2_abc(Lattice)
resolution_x = vector_a/NGX
resolution_y = vector_b/NGY
resolution_z = vector_c/NGZ
grid_pot, electrons = md.density_2_grid(vasp_pot,NGX,NGY,NGZ)

"""
Explanation: The code below is usually set in the PlanarAverage.py file; you don't need to edit it
End of explanation
"""


## POTENTIAL
planar = md.planar_average(grid_pot, NGX, NGY, NGZ)
## MACROSCOPIC AVERAGE
macro  = md.macroscopic_average(planar, lattice_vector/2, resolution_z)

"""
Explanation: The code below will prompt you to say which axis you want to average along
End of explanation
"""


fig, ax1 = plt.subplots(1, 1, sharex=True)

textsize = 22
mpl.rcParams['xtick.labelsize'] = textsize
mpl.rcParams['ytick.labelsize'] = textsize
mpl.rcParams['figure.figsize'] = (15, 8)

ax1.plot(planar,label="Planar",lw=3)
ax1.plot(macro,label="Macroscopic",lw=3)

ax1.set_ylabel('V/V', fontsize=22)
ax1.set_xlabel('Grid position', fontsize=22)

ax1.set_xlim(0,len(planar))

ax1.set_facecolor((0.95,0.95,0.95))
ax1.grid(True)

ax1.legend(fontsize=22)
plt.savefig('hj-offset.png')
plt.show()
np.savetxt(output_file,macro)

"""
Explanation: Now we can plot the results
End of explanation
"""


extrema = md.vasp_tools.get_band_extrema('OUTCAR_ZnO_av')
print(extrema)
extrema = md.vasp_tools.get_band_extrema('OUTCAR_ZnS_av')
print(extrema)

"""
Explanation: Get $\Delta V$
From inspection of the macroscopic average, saved in planar.dat we can find that the value of $V$ at each plateau is 2.64 V and -2.61 V. So the step $\Delta V$ is 5.25 V
The offset
Applying the equation from earlier, the VBO is 
1.07 - 2.28 + 5.25 = 4.04 eV
Procedure (ii)
In this procedure we can take into account the fact that the band edge positions are affected by changes to the cell volume. To do this we first calculate the eigenvalues $\epsilon_{VBM}$ at an average volume of the two materials.
End of explanation
"""


input_file = 'LOCPOT'
lattice_vector = 3.6
output_file = 'planar.dat'

vasp_pot, NGX, NGY, NGZ, Lattice = md.read_vasp_density(input_file)
vector_a,vector_b,vector_c,av,bv,cv = md.matrix_2_abc(Lattice)
resolution_x = vector_a/NGX
resolution_y = vector_b/NGY
resolution_z = vector_c/NGZ
grid_pot, electrons = md.density_2_grid(vasp_pot,NGX,NGY,NGZ)

## POTENTIAL
planar = md.planar_average(grid_pot,NGX,NGY,NGZ)
## MACROSCOPIC AVERAGE
macro  = md.macroscopic_average(planar,lattice_vector/2,resolution_z)

fig, ax1 = plt.subplots(1, 1, sharex=True)

textsize = 22
mpl.rcParams['xtick.labelsize'] = textsize
mpl.rcParams['ytick.labelsize'] = textsize
mpl.rcParams['figure.figsize'] = (10, 6)

ax1.plot(planar,label="Planar",lw=3)
ax1.plot(macro,label="Macroscopic",lw=3)

ax1.set_xlim(0,len(planar))

ax1.set_facecolor((0.95,0.95,0.95))
ax1.grid(True)

ax1.legend(fontsize=22)
plt.show()
np.savetxt(output_file,macro)

"""
Explanation: The deformation potential
To account for the effect of strain on the band positions we need to know the deformation potential $\alpha$. In this case the values for $\alpha$ of ZnO and ZnS (and many other systems) are available Phys. Rev. B 73, 245206.
$\alpha_{ZnO} = 0.48$
$\alpha_{ZnS} = 0.83$
The equation we now use is 
$ VBO = (\epsilon_{vbm}^a + \alpha_a \partial \ln V_a )- (\epsilon_{vbm}^b + \alpha_b \partial \ln V_b) + \Delta V $
Here the eigenvalues and offset potential are the same as in procedure (i), but calculated for an average volume cell. We also use the log of the change in volume of each phase multiplied by the deformation potential ($ \alpha_B \partial \ln V_B $).
The potential alignment
This procedure is the same as for procedure (i), but always using the average volume cell. Actually as we used the average cell in procedure (i), it turns out to be exactly the same. But we can do it again to reinforce.
End of explanation
"""


dlnVa = (47.55 - 61.04) / 61.04 # ZnO
dlnVb = (76.88 - 61.04) / 61.04 # ZnS

VBO_natural = (-2.464 + 0.48 * dlnVa ) - (4.8288 + 0.83 * dlnVb) + 5.25

print 'Natural offset: %3.1f eV' % (VBO_natural)

"""
Explanation: As before the value of $\Delta V$ is 5.25 V.
Volume change
By comparing the average to the equilibrium POSCARs we can get the values of $\partial\ln V = \frac{\Delta V}{V_{av}}$
End of explanation
"""


