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
#The following function reads the charge density from a cube file 
def read(cubefile):
  density = CHD()
  f = open(cubefile,'r')
  #skip two header lines
  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):

  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)
  return density

if __name__ == "__main__":
  if len(sys.argv) != 2:
    print("Incorrect number of arguments, run as ./ChargeDensity.py CUBEFILELOCATION")
  density = read(sys.argv[1])  
  #For the main function I care about the total number of electrons
# 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
    $\begingroup$ +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? $\endgroup$ Dec 4, 2022 at 21:20
  • $\begingroup$ 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. $\endgroup$
    – Sak
    Dec 4, 2022 at 21:40
  • 1
    $\begingroup$ 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. $\endgroup$ Dec 4, 2022 at 21:51
  • $\begingroup$ Sure. Thanks. @NikeDattani. And, also if you have time, please feel free to check the question if any solution idea can be generated. $\endgroup$
    – Sak
    Dec 5, 2022 at 5:43

1 Answer 1


After visualising the cube file in VMD, it was found that the coordinates of one of Si atoms (-3.9375, 3.9375, 3.9375) is beyond the 30x30x30 voxels system (interpreting 0.075 atomic units as width of a voxel from the cube file).

enter image description here

Hence, Spherical Integration is not possible as volumetric data provided seems to be insufficient. It would rather be appropriate to use Bader Charge Analysis to find total charge associated with each atom, and the zero flux surfaces which might give a better picture about the interstitial charge of the 2 Si atom system.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .