<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Ionisation potential of a porous material
In this example we use MacroDensity with VASP to align the energy levels of a porous material.
The procedure involves one DFT calculaion, yielding different important values
$\epsilon_{vbm}$ - the valence band maximum
$V_{vac}$ - the vacuum potential
The ionisation potential ($IP$) is then obtained from
Step1:  Read the potential
Step2:  Look for pore centre points
For this we will use VESTA.
Open the LOCPOT in VESTA.
Expand to 2x2x2 cell, by choosing the Boundary option on the left hand side.
Look for a pore centre - I think that [1,1,1] is at a pore centre here.
Now draw a lattice plane through that point.
Choose Edit > Lattice Planes.
Click New.
Put in the Miller index (I choose 1,1,1).
Move the plane up and down using the d parameter, until it passes through the point you think is the centre.
It should look like the picture below.
<img src="111.png">
Now we look at a contour plot of this plane to see if we are at a plateau.
Utiltiy > 2D Data Display.
Click Slice and enter the same parameters as in the 3D view.
Now choose contours to play with the settings
Z(max) and Z(min) tell you the potentail max and min.
Set contour max = Z(max) and contour min = 0
Set the interval to 0.1
With some playing with the General settings, you can get something like this
Step3:  We want to try a range of sampling area sizes.
We analyse how the potential is affects.
We also want low variance (plateau condidtions).
Ideally we should have as large an area as possible, with low (< 1e-5) variance.
Step4:  From the OUTCAR the VBM is at -2.4396 V
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
%load_ext autoreload
%autoreload 2
import sys
sys.path.append('../../')
import imp
import macrodensity as md
import math
import numpy as np
import matplotlib.pyplot as plt
import os
Explanation: Ionisation potential of a porous material
In this example we use MacroDensity with VASP to align the energy levels of a porous material.
The procedure involves one DFT calculaion, yielding different important values
$\epsilon_{vbm}$ - the valence band maximum
$V_{vac}$ - the vacuum potential
The ionisation potential ($IP$) is then obtained from:
$IP = V_{vac} - \epsilon_{vbm}$
The difference to a bulk calculation is that here the material itself has a vacuum within. That means that we can sample the vacuum level from there.
The procedure was first outlined in a seminal JACS paper, read it here.
Our system
The beautiful ZIF-8 is our porous system of choice for this demonstration.
<img src="zif8.png">
Procedure
Optimise the structre
Calculate the electronic structure at your chosen level of theory Remember in your INCAR:
LVHAR = .TRUE.  # This generates a LOCPOT file with the potential 
Locate the centre of the largest pore - do this "by eye" first
Plot the potential in that plane, so see if it plateaus
Plot a profile of the potential across the pore, again to see the plateau
Sample the potential from the pore centre
NB This whole procedure is probably better run in a notebook than by script, the reason being that you can read the file once, then do the manipulations later. The reading is the intensive and time consuming step.
End of explanation
if os.path.isfile('LOCPOT'):
    print('LOCPOT already exists')
else:
    os.system('bunzip2 LOCPOT.bz2')
input_file = 'LOCPOT'
#=== No need to edit below
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: Read the potential
End of explanation
cube_origin = [1,1,1]
travelled = [0,0,0]
int(cube_origin[0]*NGX)
Explanation: Look for pore centre points
For this we will use VESTA.
Open the LOCPOT in VESTA.
Expand to 2x2x2 cell, by choosing the Boundary option on the left hand side.
Look for a pore centre - I think that [1,1,1] is at a pore centre here.
Now draw a lattice plane through that point.
Choose Edit > Lattice Planes.
Click New.
Put in the Miller index (I choose 1,1,1).
Move the plane up and down using the d parameter, until it passes through the point you think is the centre.
It should look like the picture below.
<img src="111.png">
Now we look at a contour plot of this plane to see if we are at a plateau.
Utiltiy > 2D Data Display.
Click Slice and enter the same parameters as in the 3D view.
Now choose contours to play with the settings
Z(max) and Z(min) tell you the potentail max and min.
Set contour max = Z(max) and contour min = 0
Set the interval to 0.1
With some playing with the General settings, you can get something like this:
<img src='plane.png'>
We can see the [1,1,1], at the centre of the picture is a maximum and is a plateau, so we can now use it for sampling.
Sampling the potential
We now set the point to sample at [1,1,1]
We must also set the travelled parameter, for this type of analysis it is always [0,0,0].
End of explanation
dim = [1,10,20,40,60,80,100]
print("Dimension   Potential   Variance")
print("--------------------------------")
for d in dim:
    cube = [d, d, d]
    cube_potential, cube_var = md.volume_average(cube_origin, cube,grid_pot, NGX, NGY, NGZ, travelled=travelled)
    print(" %3i     %10.4f   %10.6f"%(d, cube_potential, cube_var))
Explanation: We want to try a range of sampling area sizes.
We analyse how the potential is affects.
We also want low variance (plateau condidtions).
Ideally we should have as large an area as possible, with low (< 1e-5) variance.
End of explanation
print("IP: %3.4f eV" % (2.3068 -- 2.4396 ))
Explanation: From the OUTCAR the VBM is at -2.4396 V
End of explanation
<END_TASK>