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() 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()) #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) 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*density.N*density.N)) #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) #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.