Coverage for gpaw/analyse/multipole.py: 98%

1import numpy as np 

2 

3from ase.units import Bohr 

4from ase.parallel import paropen 

5 

6from gpaw.spherical_harmonics import Y 

7from gpaw.utilities.tools import coordinates 

8 

9 

10class Multipole: 

11 

12 """Expand a function on the grid in multipole moments 

13 relative to a given center. 

14 

15 center: Vector [Angstrom] 

16 """ 

17 

18 def __init__(self, center, calculator=None, lmax=6): 

19 

20 self.center = center / Bohr 

21 self.lmax = lmax 

22 

23 self.gd = None 

24 self.y_Lg = None 

25 self.l_L = None 

26 

27 if calculator is not None: 

28 self.initialize(calculator.density.finegd) 

29 

30 def initialize(self, gd): 

31 """Initialize Y_L arrays""" 

32 

33 self.gd = gd 

34 

35 r_cg, r2_g = coordinates(gd, self.center, tiny=1.e-78) 

36 r_g = np.sqrt(r2_g) 

37 rhat_cg = r_cg / r_g 

38 

39 self.l_L = [] 

40 self.y_Lg = [] 

41 npY = np.vectorize(Y, (float,), 'spherical harmonic') 

42 L = 0 

43 for l in range(self.lmax + 1): 

44 for m in range(2 * l + 1): 

45 self.y_Lg.append( 

46 np.sqrt(4 * np.pi / (2 * l + 1)) * r_g ** l * 

47 npY(L, rhat_cg[0], rhat_cg[1], rhat_cg[2]) 

48 ) 

49 self.l_L.append(l) 

50 L += 1 

51 

52 def expand(self, f_g): 

53 """Expand a function f_g in multipole moments 

54 units [e * Angstrom**l]""" 

55 

56 assert f_g.shape == self.gd.empty().shape 

57 

58 q_L = [] 

59 for L, y_g in enumerate(self.y_Lg): 

60 q_L.append(self.gd.integrate(f_g * y_g)) 

61 q_L[L] *= Bohr ** self.l_L[L] 

62 

63 return np.array(q_L) 

64 

65 def to_file(self, calculator, 

66 filename='multipole.dat', 

67 mode='a'): 

68 """Expand the charge distribution in multipoles and write 

69 the result to a file""" 

70 

71 if self.gd is None: 

72 self.initialize(calculator.density.finegd) 

73 q_L = self.expand(-calculator.density.rhot_g) 

74 

75 f = paropen(filename, mode) 

76 

77 print('# Multipole expansion of the charge density', file=f) 

78 print('# center =', self.center * Bohr, 'Angstrom', file=f) 

79 print('# lmax =', self.lmax, file=f) 

80 print('# see https://gitlab.com/gpaw/gpaw/-/blob/master/c/bmgs/' 

81 'sharmonic.py', file=f) 

82 print('# for the definition of spherical harmonics', file=f) 

83 print('# l m q_lm[|e| Angstrom**l]', file=f) 

84 

85 L = 0 

86 for l in range(self.lmax + 1): 

87 for m in range(-l, l + 1): 

88 print(f'{l:2d} {m:3d} {q_L[L]:g}', file=f) 

89 L += 1 

90 f.close()