Tag: Bader

GPAW installation with pip

Between installation with conda and compilation of libraries, an intermediate path – installation of GPAW with pip – is a compromise for those who wish to text specific GPAW branches or packages.

For example, I wish to text self-interaction error correction (SIC) and evaluate Bader charges with pybader. Neither SIC nor pybader is compatible with the recent GPAW. Here is not to get a workable version.

# numba in pybader is not compatible with python 3.11, so create a conda environment with python 3.10
conda create -n gpaw-pip python=3.10 
conda activate gpaw-pip

conda install -c conda-forge libxc libvdwxc
conda install -c conda-forge ase
# ensure that you install the right openmpi (not external)
conda install -c conda-forge openmpi ucx
conda install -c conda-forge compilers
conda install -c conda-forge openblas scalapack
conda install -c conda-forge pytest
pip install pybader

# Get a developer version of GPAW with SIC
git clone -b dm_sic_mom_update https://gitlab.com/alxvov/gpaw.git
cd gpaw
cp siteconfig_example.py siteconfig.py

# In the siteconfig.py rewrite
'''
fftw = True
scalapack = True
if scalapack:
    libraries += ['scalapack']
'''

unset CC
python -m pip install -e .
gpaw info

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