I have a cube file of charge density. The unit cell consists two silicon atoms.

If I use a python code to calculate the net charge of the system, it results around 24, which is ok because the pseudopotential of each silicon atom consists 12 electrons. The code is given at the end, if needed to see.

Now, I would request your help to have a python code that can actually calculate charge of an individual Si atom from that cube file (rather than considering both the atoms)? I think the charge should be around 12 for this case (for individual atom). May be a spherical integration up to a certain cutoff limit can help (I am not sure though about how to select that limit/radius). Please help, if you have some time.

Here is the cube file: 3D charge density

Here is the python code that can calculate the total charge from that cube file:

#!/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 runs as

python ChargeDensity.py si_pseudo.cube

2 Answers 2


If you have a .CUB or .CUBE file already then it is as simple as using this python library- CubeToolz

It can perform the following operations:

  • Read and write Gaussian cube files
  • Translate and rotate cube data
  • Integrate around a particular atom
  • Integrate around a sphere
  • Integrate around the whole cube file
  • Take the planar average Add/subtract/multiply cube files

The relevant option for you is to integrate the whole cube file.

cube_tools.py -i file.cube

This is one of my go-to libraries for analyzing densities in .CUBE format.

  • $\begingroup$ Thanks a lot. I will surely check it. @Newbie $\endgroup$
    – Sak
    Dec 9, 2022 at 19:06

I do not have a python code for you. But I could direct 2 post-processing open software you could use to calculate individual charge of your atoms.

(1) MultiWfn (http://sobereva.com/multiwfn/)

(2) Bader Charge Analysis (http://theory.cm.utexas.edu/henkelman/code/bader/)

Bader Charge Analysis is based on charge density partitioning developed by, you guessed it, Bader. Multiwfn can do other charge integration based on different partitioning such as tessellations.

I suggest trying them out with tutorial samples before applying it for your system.


  • $\begingroup$ The user is asking to add charge density upto custom radius of sphere region in space. An atom is situated at some point (x, y, z) what will be total charge in a spherical region of radius (r). Bader divides charge based on gradient to individual atoms one need not define radius. $\endgroup$ Dec 9, 2022 at 11:59
  • $\begingroup$ @Pranavkumar As I read it as the user isn't clear what they want - "May be a spherical integration up to a certain cutoff limit can help (I am not sure though about how to select that limit/radius)" (emphasis mine). As such pointing out established methods such as Bader analysis to perform what they require could be extremely useful. $\endgroup$
    – Ian Bush
    Dec 9, 2022 at 14:59
  • $\begingroup$ Thank you very much @Pranavkumar, yes you are right $\endgroup$
    – Sak
    Dec 9, 2022 at 19:04
  • $\begingroup$ Thanks Ian Bush, yes Bader analyses might help. The current problem (goal) would be more meaningful with python code though. $\endgroup$
    – Sak
    Dec 9, 2022 at 19:06

You must log in to answer this question.

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