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

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
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()[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)
#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

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.

• Thanks a lot. I will surely check it. @Newbie
– Sak
Commented 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/)

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!

• 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. Commented Dec 9, 2022 at 11:59
• @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. Commented Dec 9, 2022 at 14:59
• Thank you very much @Pranavkumar, yes you are right
– Sak
Commented Dec 9, 2022 at 19:04
• Thanks Ian Bush, yes Bader analyses might help. The current problem (goal) would be more meaningful with python code though.
– Sak
Commented Dec 9, 2022 at 19:06