Bader and Hirshfeld charges with python

Bader analysis is a fast and simple way of getting atomic charges. It is especially useful for periodic calculations. The analysis can be done on the fly using pybader tool from pypi.org/project/pybader. I recommend using it within conda environments.

The installation is straightforward:

pip install pybader

The usage is less obvious. Here is a function for obtaining xmol-type xyz that can obtained with GPAW and visualized with ASE:


def xyzb(atoms, filename, nCPU):
  from pybader.io import gpaw
  from pybader.interface import Bader
  import os
  bader = Bader(*gpaw.read_obj(atoms.calc)) # read ASE object 'atoms'
  bader(threads=nCPU)                       # specify the number of CPUs
  f = open('{0}.xyz'.format(filename), 'w') # set xmol format
  b = bader.atoms_charge                    # get number of electrons per atom
  n = atoms.get_atomic_numbers()            # get atomic numbers
  a = atoms.get_chemical_symbols()          # get chemical symbols
  p = atoms.get_positions()                 # get positions of the atoms
  f.write('{0}\n'.format(len(a)))
  f.write('Properties=species:S:1:pos:R:3:charge:R:1\n') # ensure compatibility with ASE
  for i in range(len(a)):                   # print symbol, positions, and charge
    s = '{0}'.format(a[i])
    x = '{0:.6f}'.format(round(p[i][0],6))
    y = '{0:.6f}'.format(round(p[i][1],6))
    z = '{0:.6f}'.format(round(p[i][2],6))
    c = '{0:.3f}'.format(round(float(n[i]) - float(b[i]),3))
    f.write('{0:<4}{1:>16}{2:>16}{3:>16}{4:>10}\n'.format(s,x,y,z,c))
  f.close()
  del bader
  os.remove('bader.p')

Similarly one can obtain xmol-type xyz file with Hirshfeld charges:


def xyzh(atoms, filename):
  from gpaw.analyse.hirshfeld import HirshfeldPartitioning
  f = open('{0}.xyz'.format(filename), 'w') # set xmol format
  a = atoms.get_chemical_symbols()          # get chemical symbols
  p = atoms.get_positions()                 # get positions of the atoms
  h = HirshfeldPartitioning(atoms.calc).get_charges()
  f.write('{0}\n'.format(len(a)))
  f.write('Properties=species:S:1:pos:R:3:charge:R:1\n') # ensure compatibility with ASE
  for i in range(len(a)):                   # print symbol, positions, and charge
    s = '{0}'.format(a[i])
    x = '{0:.6f}'.format(round(p[i][0],6))
    y = '{0:.6f}'.format(round(p[i][1],6))
    z = '{0:.6f}'.format(round(p[i][2],6))
    c = '{0:.3f}'.format(round(h[i],3))
    f.write('{0:<4}{1:>16}{2:>16}{3:>16}{4:>10}\n'.format(s,x,y,z,c))
  f.close()

In the LCAO mode of GPAW one can also get the Mulliken charges. Test before using:


def xyzm(atoms, filename):
  from gpaw.lcao.tools import get_mulliken
  f = open('{0}.xyz'.format(filename), 'w') # set xmol format
  a = atoms.get_chemical_symbols()          # get chemical symbols
  p = atoms.get_positions()                 # get positions of the atoms
  m = get_mulliken(atoms.calc, range(len(a)))
  f.write('{0}\n'.format(len(a)))
  f.write('Properties=species:S:1:pos:R:3:charge:R:1\n') # ensure compatibility with ASE
  for i in range(len(a)):                   # print symbol, positions, and charge
    s = '{0}'.format(a[i])
    x = '{0:.6f}'.format(round(p[i][0],6))
    y = '{0:.6f}'.format(round(p[i][1],6))
    z = '{0:.6f}'.format(round(p[i][2],6))
    c = '{0:.3f}'.format(round(m[i],3))
    f.write('{0:<4}{1:>16}{2:>16}{3:>16}{4:>10}\n'.format(s,x,y,z,c))
  f.close()

ILMAT5 presentation

LINK to the PRESENTATION

Abstract

In this work [0], we have applied the DFT-based delta Kohn–Sham (ΔKS) method to ion pairs in a vacuum to obtain X-ray pho­toelectron spectra of corresponding ionic liquids (IL). On the example of forty ion pairs, we demonstrate how the core level binding energy (BE) values can be calcu­lated and used to plot theo­retical spectra at a low computational cost. Furthermore, we compare the ΔKS results, 1s Kohn–Sham orbital energies, and atomic charges against the experi­mental X-ray photoelec­tron data. Recently, in connection to the electro­chemical application in the super­capacitors, we have measured spectra for EMImBF4 and EMImB(CN)4 ionic liquids at the carbon–IL interface [1–3]. Other experimental spectra were obtained from the literature [4,5]. Both the ΔKS BE values and the 1s Kohn–Sham orbital energies show a correlation, yet with a different order of the BEs assigned to spe­cific atoms. We find that neither DDEC6 nor Bader charges cor­relate with the experi­mental data. Thus, the DFT calculations of 1s Kohn–Sham orbital energies provide the fastest way of pre­dicting the XPS spectra. However, more detailed experimental studies are required to resolve the right order of the BE values and its rela­tion to the atomistic structure of the ILs. The ΔKS calculations provide the most precise estimations of the BEs. Herewith, they also demand more resources and cause computa­tional difficulties discussed in the presenta­tion. Besides the prediction power, a robust computational method can help in intepre­tating experimental data when the appropriate reference values are either not available nor directly applicable. Thus, the ΔKS method can find its application in various fields of physics and chemistry where the XPS is used for re­solving electronic and geometric structures of pure ILs, their mixtures, and at interfaces.

In this work, we have applied the DFT-based delta Kohn–Sham (ΔKS) method to ion pairs in a vacuum to obtain X-ray pho­toelectron spectra of corresponding ionic liquids (IL). On the example of forty ion pairs, we demonstrate how the core level binding energy (BE) values can be calcu­lated and used to plot theo­retical spectra at a low computational cost. Furthermore, we compare the ΔKS results, 1s Kohn–Sham orbital energies, and atomic charges against the experi­mental X-ray photoelec­tron data. Recently, in connection to the electro­chemical application in the super­capacitors, we have measured spectra for EMImBF4 and EMImB(CN)4 ionic liquids at the carbon–IL interface [1–3]. Other experimental spectra were obtained from the literature [4,5]. Both the ΔKS BE values and the 1s Kohn–Sham orbital energies show a correlation, yet with a different order of the BEs assigned to spe­cific atoms. We find that neither DDEC6 nor Bader charges cor­relate with the experi­mental data. Thus, the DFT calculations of 1s Kohn–Sham orbital energies provide the fastest way of pre­dicting the XPS spectra. However, more detailed experimental studies are required to resolve the right order of the BE values and its rela­tion to the atomistic structure of the ILs. The ΔKS calculations provide the most precise estimations of the BEs. Herewith, they also demand more resources and cause computa­tional difficulties discussed in the presenta­tion. Besides the prediction power, a robust computational method can help in intepre­tating experimental data when the appropriate reference values are either not available nor directly applicable. Thus, the ΔKS method can find its application in various fields of physics and chemistry where the XPS is used for re­solving electronic and geometric structures of pure ILs, their mixtures, and at interfaces.

[0] M. Lembinen, E. Nõmmiste, H. Ers, B. Docampo‐Álvarez, J. Kruusma, E. Lust, V.B. Ivaništšev, Calculation of core‐level electron spectra of ionic liquids, Int. J. Quantum Chem. 120 (2020). https://doi.org/10.1002/qua.26247.

[1] J. Kruusma, A. Tõnisoo, R. Pärna, E. Nõmmiste, I. Tallo, T. Romann, E. Lust, Electrochimica Acta 206 (2016) 419–426.

[2] J. Kruusma, A. Tõnisoo, R. Pärna, E. Nõmmiste, I. Kuusik, M. Vahtrus, I. Tallo, T. Romann, E. Lust, J. Electrochem. Soc. 164 (2017) A3393–A3402.

[3] A. Tõnisoo, J. Kruusma, R. Pärna, A. Kikas, M. Hirsimäki, E. Nõmmiste, E. Lust, J. Electrochem. Soc. 160 (2013) A1084–A1093.

[4] A. Foelske-Schmitz, D. Weingarth, R. Kötz, Surf. Sci. 605 (2011) 1979–1985.

[5] I.J. Villar-Garcia, E.F. Smith, A.W. Taylor, F. Qiu, K.R.J. Lovelock, R.G. Jones, P. Licence, Phys. Chem. Chem. Phys. 13 (2011) 2797–2808.