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()
