5
$\begingroup$

I have calculated the charge densities (in 1D, 2D and 3D formats) with the help of Quantum Espresso DFT.

The system (/material) is simple, two atoms in an unit cell.

However, I want to determine the number of electrons (in the interstitial position between the two atoms). Is there any python code with numerical approximation methods anyone can guide? It would be great to get your detailed suggestion.

For example, the following image is for the charge densities (in atomic units) of that material under different conditions. I need to know the number of electrons near the region 2 alat (where there is interstitial bumps), which is in between two atoms (two black dotted vertical lines), one reference atom and it's nearest neighbor atom. The image is plotted from an 1D data set of charge density vs distance for different structures.

This is from the 1D charge density data

More information (optional):

Lattice constant (for one structure): 4.5 Bohr

Ref. atom   x   y   z

alat    0.375   0.125   0.625

Bohr    1.6875  0.5625  2.8125

Angstrom    0.89298619  0.297662    1.48831031

Nearest neighbor atom

alat    0.625   0.875   0.375

Bohr    2.8125  3.9375  1.6875

Angstrom  1.48831031    2.083634    0.89298619

The interstitial position might be considered maybe halfway between the two atoms (I am not sure though).

Here is the 3D charge density file for same structure (for which the lattice constant is given above): 3D cube file

The total charge (number of electrons) can be calculated from here:

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

This python code runs as (and gives output as 24, which is the total charge of the unit cell):

python ChargeDensity.py si_pseudo.cube

So again, we need to modify the python code in a way so that it can determine the interstitial charge from the 3D cube file of charge density. Hope to get your cordial cooperation.

$\endgroup$
3
  • 1
    $\begingroup$ In principle, your need: "number of electrons" make no sense. As the electron is considered as a wave, it is spread over space, so, it is not localized in a sense of a particle. You can make QTAIM analysis where you can classify the critical points and have an idea of how intense is the electron sharing between the two cores. $\endgroup$
    – Camps
    Nov 9, 2022 at 23:26
  • $\begingroup$ Thanks @Camps for your response. The bader analysis is something like determining number of interstitial electrons. It can be done in different other approaches as well. The question was to know about any of those steps. $\endgroup$
    – Sak
    Nov 10, 2022 at 18:52
  • $\begingroup$ I gave my +1 long ago, but have you managed to figure this out over the last 6+ months? Are you still actively needing an answer? Please update us! $\endgroup$ May 28, 2023 at 3:42

0