1from pathlib import Path
2from typing import Union
3
4import numpy as np
5
6
7def read_bader_charges(filename: Union[str, Path] = 'ACF.dat') -> np.ndarray:
8    path = Path(filename)
9    charges = []
10    with path.open() as fd:
11        for line in fd:
12            words = line.split()
13            if len(words) == 7:
14                charges.append(float(words[4]))
15    return np.array(charges)
16
17
18if __name__ == '__main__':
19    import subprocess
20    import sys
21    from ase.io import write
22    from ase.units import Bohr
23    from gpaw import GPAW
24    from gpaw.utilities.ps2ae import PS2AE
25
26    calc = GPAW(sys.argv[1])
27    converter = PS2AE(calc, grid_spacing=0.05)
28    density = converter.get_pseudo_density()
29    ne = density.sum() * converter.dv
30    print(ne, 'electrons')
31    write('density.cube', calc.atoms, data=density * Bohr**3)
32    subprocess.run('bader -p all_atom density.cube'.split())
33    charges = read_bader_charges()
34    for setup, charge in zip(calc.density.setups, charges):
35        charge -= setup.Nv
36        print(f'{setup.symbol:2} {charge:10.6f}')
37