Hello All,
I hope you are starting to enjoy summer as much as I am.
This is a quick tutorial on how to read a cube file into a numpy array in python and integrate it. This was born out of a very unique problem I am having right now involving adding electrons to a system to see the effects on the fermi surface in wide bandgap semiconductors. I will get to the bottom of it at some point, but I figured this would be a cool little tutorial to save you some coding in the future.
The first step is to generate a cube file from quantum espresso. Quantum Espresso stores the charge density of a calculation in a binary file called charge-density.dat (that you can find in your *.save directory). You could, and I have before, read the charge density directly from that file, however binary files are unwieldy and anything that I make that can post-process that would be specific to QE. A cube file has a standard format that many different programs generate. You can see a basic outline of how this file format can be found here:http://paulbourke.net/dataformats/cube/
To generate a cube file of the total charge density, you will need to use pp.x. This utility simply does the post processing of the binary file that I was trying to avoid above. By using the following example input file, you can run: pp.x < pp.in > pp.out to generate two files: al.filplot and chd.cube. This is the total charge density of the system; the file we are interested for this tutorial is the chd.cube file.
&INPUTPP prefix='al' outdir='.' filplot='al.pp' plot_num = 0 / &PLOT nfile=1 filepp(1)='al.pp' iflag=3 output_format=6 fileout='chd.cube' /
Now, say we are interested in reading the chd.cube file and integrating to determine the total number of electrons in our system, how would we do it? Well, for starters, the charge density is stored at discrete volumes of dV in a grid of dimensions of Nx x Ny x Nz. If we read the charge density into a numpy array and sum up the array we can determine the number of electrons by multiplying this by the volume element dV (discrete integration). I have written this up into another python module (I am turning into quite the python hipster), that you can directly download here: http://www.levilentz.com/Codes/ChargeDensity.py
I have also pasted it here for posterity. Essentially by running it as ./ChargeDensity.py chd.cube it will read the file into memory and perform the integration. Looking at my problem, I tested it on a positively charged system, and sure enough there is an extra half an electron in the system. Which is leaving me scratching my head more than anything; I guess that is just research.
Happy computing,
Levi
#!/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([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())