How to calculate the individual atomic charge from a cube file of charge density?

Question

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
    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 runs as

python ChargeDensity.py si_pseudo.cube

Answer 1

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.

Answer 2

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.

Cheers!