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.

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.



In [1]:

    
%matplotlib inline
%load_ext autoreload
%autoreload 2
import imp
import macrodensity as md
import math
import numpy as np
import matplotlib.pyplot as plt

Read the potential



In [2]:

    
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)









    



Reading potential at point 100000
Reading potential at point 200000
Reading potential at point 300000
Reading potential at point 400000
Reading potential at point 500000
Reading potential at point 600000
Reading potential at point 700000
Reading potential at point 800000
Reading potential at point 900000
Reading potential at point 1000000
Reading potential at point 1100000
Reading potential at point 1200000
Reading potential at point 1300000
Reading potential at point 1400000
Reading potential at point 1500000
Reading potential at point 1600000
Reading potential at point 1700000
Reading potential at point 1800000
Reading potential at point 1900000
Reading potential at point 2000000
Reading potential at point 2100000
Reading potential at point 2200000
Reading potential at point 2300000
Reading potential at point 2400000
Reading potential at point 2500000
Reading potential at point 2600000
Reading potential at point 2700000
Reading potential at point 2800000
Reading potential at point 2900000
Reading potential at point 3000000
Reading potential at point 3100000
Reading potential at point 3200000
Reading potential at point 3300000
Reading potential at point 3400000
Reading potential at point 3500000
Reading potential at point 3600000
Reading potential at point 3700000
Reading potential at point 3800000
Reading potential at point 3900000
Reading potential at point 4000000
Reading potential at point 4100000
Reading potential at point 4200000
Reading potential at point 4300000
Reading potential at point 4400000
Reading potential at point 4500000
Reading potential at point 4600000
Reading potential at point 4700000
Reading potential at point 4800000
Reading potential at point 4900000
Reading potential at point 5000000
Reading potential at point 5100000
Reading potential at point 5200000
Reading potential at point 5300000
Reading potential at point 5400000
Reading potential at point 5500000
Reading potential at point 5600000
Reading potential at point 5700000
Reading potential at point 5800000
Reading potential at point 5900000
Reading potential at point 6000000
Reading potential at point 6100000
Reading potential at point 6200000
Reading potential at point 6300000
Reading potential at point 6400000
Reading potential at point 6500000
Reading potential at point 6600000
Reading potential at point 6700000
Reading potential at point 6800000
Reading potential at point 6900000
Reading potential at point 7000000
Reading potential at point 7100000
Reading potential at point 7200000
Reading potential at point 7300000
Reading potential at point 7400000
Reading potential at point 7500000
Reading potential at point 7600000
Reading potential at point 7700000
Reading potential at point 7800000
Reading potential at point 7900000
Reading potential at point 8000000
Reading potential at point 8100000
Reading potential at point 8200000
Reading potential at point 8300000
Reading potential at point 8400000
Reading potential at point 8500000
Reading potential at point 8600000
Reading potential at point 8700000
Reading potential at point 8800000
Reading potential at point 8900000
Reading potential at point 9000000
Reading potential at point 9100000
Reading potential at point 9200000
Reading potential at point 9300000
Reading potential at point 9400000
Reading potential at point 9500000
Reading potential at point 9600000
Reading potential at point 9700000
Reading potential at point 9800000
Reading potential at point 9900000
Reading potential at point 10000000
BBBB		OOOO		OOOO		MMMMM	
BBBB		OOOO		OOOO		MMMMM	
BBBB		OOOO		OOOO		MMMMM	
B  B	        OOOO		OOOO		MMMMM	
B  B	        O  O		O  O		MMMMM	
B  B	        O  O		O  O		MMMMM	
B  B	        O  O		O  O		MMMMM	
B  B	        O  O		O  O		MMMMM	
BBBB	        O  O            O  O		M M M	
BBBB	        O  O		O  O		M M M	
BBBB	        O  O		O  O		M M M	
B  B	        O  O		O  O		M M M	
B  B	        O  O		O  O		M M M	
B  B	        O  O		O  O		M M M	
B  B	        O  O		O  O		M M M	
B  B	        OOOO    	OOOO		M M M	
BBBB            OOOO	        OOOO		M M M	
BBBB            OOOO	        OOOO		M M M	
BBBB            OOOO	        OOOO		M M M	
('Average of the potential = ', -2.2891856197896633e-07)

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.
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:

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].



In [3]:

    
cube_origin = [1,1,1]
travelled = [0,0,0]

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.



In [4]:

    
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.cube_potential(cube_origin,travelled,cube,grid_pot,NGX,NGY,NGZ)
    print " %3i     %10.4f   %10.6f"%(d,cube_potential,cube_var)









    



Dimension   Potential   Variance
--------------------------------
   1         2.3068     0.000000
  10         2.3068     0.000001
  20         2.3068     0.000003
  40         2.3068     0.000019
  60         2.3067     0.000108
  80         2.3048     0.001151
 100         2.2883     0.015872



In [5]:

    
extrema = md.vasp_tools.get_band_extrema('OUTCAR')
print extrema









    



[-2.4396, 3.4777]



In [6]:

    
print "IP: %3.4f eV" % (2.3068 -- 2.4396 )









    



IP: 4.7464 eV

AFI

No wtry the procedure with the files in the AFI folder



In [ ]: