# Interstitial charge from a cube file: Need help in python coding

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:

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
density = CHD()
f = open(cubefile,'r')
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)
#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.

• +1. I know I mentioned this on one of your previous questions, but the same is true for this question. I can't open that file from the UC Merced BOX folder. Are you able to create a folder called 10048 (the number in the URL of this question) in this repository? If not can you let me know so that I can help you? Dec 4, 2022 at 21:20
• Hi @NikeDattani, Thanks for your suggestion. Somehow could not create that folder. I tried to upload the file in that repository still. Could you please check if it's ok or please show me how to create directory in there to upload file. Thanks again.
– Sak
Dec 4, 2022 at 21:40
• Okay I fixed it. It would be helpful if you were to find a tutorial on Git and to work through the exercises to learn how to do some basics, like adding files to folders. Dec 4, 2022 at 21:51
• Sure. Thanks. @NikeDattani. And, also if you have time, please feel free to check the question if any solution idea can be generated.
– Sak
Dec 5, 2022 at 5:43 