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