#!/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())