Skip to content

Commit

Permalink
Script to sum up the charge density around an atom
Browse files Browse the repository at this point in the history
  • Loading branch information
Asif-Iqbal-Bhatti authored Jul 16, 2023
1 parent 2028622 commit 3d0dca9
Showing 1 changed file with 45 additions and 0 deletions.
45 changes: 45 additions & 0 deletions get_TotEle.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,45 @@
#!/usr/bin/env python3

'''
#=====================================================
# AUTHOR: Asif Iqbal
# GitHUB: Asif_em
# PURPOSE: Read CHGCAR file and compute the integrated
# charge density to sum up to # of electrons.
# Dependency on pymatgen module
# The following code is uploaded to my GitHub
#
#=====================================================
'''

import numpy as np
import matplotlib.pyplot as plt
from pymatgen.io.vasp.outputs import VolumetricData

poscar, data, data_aug = VolumetricData.parse_file('CHGCAR')
latt = np.array(poscar.structure.lattice.matrix)

# As per VASP the CHGCAR needs to be divided by Volume
# COuld be total or diff if ISPIN = 2
chgd = data["total"] / np.linalg.det(latt)
pos = poscar.structure.cart_coords
# get volume per grid point
dv = np.linalg.det(latt)/(np.array(chgd).shape[0]*np.array(chgd).shape[1]*np.array(chgd).shape[2])
# get the step size in each direction (x, y, z)
dr = latt/np.array(chgd).shape
total_elect = np.sum(np.array(chgd))*dv
num_grid_points = np.prod(np.array(chgd).shape)
integrated_charge = np.cumsum(np.array(chgd)) * dv
# Measure the charge density sphere
distance = np.linspace(0, np.linalg.norm(dr), num_grid_points)

print(dv, np.linalg.norm(dr))
print(f"Shape of data in CHGCAR file:: {np.array(chgd).shape}")
print(f"Total integrated electron number = {total_elect:8.5f}")
# ploting
plt.plot(distance, integrated_charge, '-o', color='m', linewidth=0.5, markersize=0.5, markerfacecolor='none')
plt.xlabel("Distance from the freestanding O atom (Å)")
plt.ylabel("Integrated charge")
plt.savefig('dvsQ.png', dpi=400, bbox_inches='tight')
plt.show()

0 comments on commit 3d0dca9

Please sign in to comment.