Coverage for gpaw/analyse/eed.py: 95%

1import numpy as np 

2from ase.units import Bohr, Hartree 

3from ase.parallel import paropen 

4from ase.data import vdw_radii 

5 

6import gpaw.cgpaw as cgpaw 

7from gpaw.io.fmf import FMF 

8 

9 

10class ExteriorElectronDensity: 

11 

12 """Exterior electron density to describe MIES spectra. 

13 

14 Simple approach to describe MIES spectra after 

15 Y. Harada et al., Chem. Rev. 97 (1997) 1897 

16 """ 

17 

18 def __init__(self, gd, atoms): 

19 """Find the grid points outside of the van der Waals radii 

20 of the atoms""" 

21 

22 assert gd.orthogonal 

23 self.gd = gd 

24 

25 n = len(atoms) 

26 atom_c = atoms.positions / Bohr 

27 vdWradius = np.empty(n) 

28 for a, atom in enumerate(atoms): 

29 vdWradius[a] = self.get_vdWradius(atom.number) 

30 

31 # define the exterior region mask 

32 mask = gd.empty(dtype=int) 

33 cgpaw.eed_region(mask, atom_c, gd.beg_c, gd.end_c, 

34 gd.h_cv.diagonal().copy(), vdWradius) 

35 self.mask = mask 

36 

37 def get_weight(self, psit_G): 

38 """Get the weight of a wave function in the exterior region 

39 (outside of the van der Waals radius). The augmentation sphere 

40 is assumed to be smaller than the van der Waals radius and hence 

41 does not contribute.""" 

42 

43 # smooth part 

44 weigth = self.gd.integrate(np.where(self.mask == 1, 

45 (psit_G * psit_G.conj()).real, 

46 0.0)) 

47 

48 return weigth 

49 

50 def get_vdWradius(self, Z): 

51 """Return van der Waals radius in Bohr""" 

52 r = vdw_radii[Z] / Bohr 

53 if np.isnan(r): 

54 msg = 'van der Waals radius for Z=' + str(Z) + ' not known!' 

55 raise RuntimeError(msg) 

56 else: 

57 return r 

58 

59 def write_mies_weights(self, wfs): 

60 file = 'eed_mies.dat' 

61 with paropen(file, 'a') as out: 

62 fmf = FMF(['exterior electron density weights after', 

63 'Y. Harada et al., Chem. Rev. 97 (1997) 1897']) 

64 print(fmf.header(), end=' ', file=out) 

65 print(fmf.data(['band index: n', 

66 'k-point index: k', 

67 'spin index: s', 

68 'k-point weight: weight', 

69 'energy: energy [eV]', 

70 'occupation number: occ', 

71 'relative EED weight: eed_weight']), end=' ', 

72 file=out) 

73 

74 print( 

75 '#; n k s weight energy occ eed_weight', 

76 file=out) 

77 for kpt in wfs.kpt_u: 

78 for n in range(wfs.bd.nbands): 

79 print('%4d %3d %1d %8.5f %10.5f %10.5f %10.5f' % 

80 (n, kpt.k, kpt.s, kpt.weight, 

81 kpt.eps_n[n] * Hartree, 

82 kpt.f_n[n], 

83 self.get_weight(kpt.psit_nG[n]) 

84 ), file=out) 

85 if hasattr(out, 'flush'): 

86 out.flush()