I have a test material- two silicon atoms in an unit cell. This particular structure shows evidence of having interstitial charge between the two atoms. I need to find out the amount of this interstitial charge. I have the cube file of the charge density (name: si_pseudo.cube). So far I have been able to find out the total charge (24, twelve electrons from each atom) in the system by a python script (ChargeDensity.py) with command:
Could you please answer in details about the python script to determine the interstitial charge?
The interstitial charge needs to be determined using a suitable numerical approximation method (may be spherical integration).
The following python script reads the charge density file and determines the total charge. But I need to determine the interstitial charge (it is usually found in the mid position between the two atoms).
#!/usr/bin/env python
import numpy as np
import sys
class CHD():
def __init__(self):
#This simply allocates the different data structures we need:
self.natoms = 0
self.grid = np.zeros(0)
self.v = np.zeros([3,3])
self.N = np.zeros([3])
self.dV = 0
def set_dV(self):
#The charge density is stored per volume. If we want to integrate the charge density
#we need to know the size of the differential volume
self.dV = 0
x = self.v[:,0]
y = self.v[:,1]
z = self.v[:,2]
self.dV = np.dot(x,np.cross(y,z))
def integrate(self):
#This allows us to integrate the stored charge density
return(np.sum(self.grid)*self.dV)
#The following function reads the charge density from a cube file
def read(cubefile):
density = CHD()
f = open(cubefile,'r')
#skip two header lines
next(f)
next(f)
line = next(f)
#Get the number of atoms if we want to store it
density.natoms = int(line.split()[0])
#This gets the nx,ny,nz info of the charge density
#As well as the differential volume
for i in range(0,3):
line = next(f).split()
density.N[i] = int(line[0])
for j in range(1,4):
density.v[i][j-1] = float(line[j])
#As of now we dont care about the positions of the atoms,
#But if you did you could read them here:
for i in range(0,density.natoms):
next(f)
density.set_dV()
density.grid = np.zeros(int(density.N[0]*density.N[1]*density.N[2]))
#This reads the data into a 1D array of size nx*ny*nz
count = 0
for i in f:
for j in i.split():
density.grid[count] = float(j)
count+=1
f.close()
return density
if __name__ == "__main__":
if len(sys.argv) != 2:
print("Incorrect number of arguments, run as ./ChargeDensity.py CUBEFILELOCATION")
sys.exit(6)
density = read(sys.argv[1])
#For the main function I care about the total number of electrons
print(density.integrate())
# Code source (from Dr. Levi): https://github.com/levilentz
The above python script worked for me as:
python ChargeDensity.py si_pseudo.cube
And it results a value: 24, which is the total charge of the unit cell)
The charge density cube file is here.
