Coverage for gpaw/utilities/ewald.py: 97%

1from math import pi 

2 

3import numpy as np 

4from scipy.special import erf, erfc 

5 

6 

7class Ewald: 

8 """Class for calculating Ewald summations. 

9 

10 'cell' is the unit cell in the cartesian basis. 

11 

12 'G' is a normalized Ewald parameter. Large G results in fast 

13 real-space convergence but slow convergence in reciprocal space. 

14 This describes the width (in reciprocal space, and inverse width 

15 in real space) of the Gaussian shaped probe charge. 

16 

17 Ng and Nl are lists specifying the number or nearest neighbors in 

18 sums over the reciprocal lattice and the real space lattice 

19 respectively. 

20 """ 

21 

22 def __init__(self, cell, G=5, Ng=[9, 9, 9], Nl=[3, 3, 3]): 

23 self.cell = cell 

24 self.Vcell = abs(np.linalg.det(cell)) 

25 self.recip_cell = np.linalg.inv(self.cell).T 

26 self.Ng = Ng 

27 self.Nl = Nl 

28 self.G = G 

29 

30 # Renormalize G to the longest lattice vector 

31 self.G /= np.sqrt((cell**2).sum(1).max()) 

32 

33 def get_wigner_seitz_radius(self, N): 

34 """Wigner-Seitz radius for N electrons in the unit cell.""" 

35 return (3 * self.Vcell / (4 * pi * N))**(1 / 3) 

36 

37 def get_sum_recip_ij(self, r_v, eps=1e-10): 

38 r"""The reciprocal space sum. 

39 

40 :: 

41 

42 ----- -x i g.r 

43 pi \ e e 

44 ---2 | ------- , with x = g^2/(4 G^2) 

45 V G / x 

46 ----- 

47 g not 0 

48 """ 

49 N = self.Ng 

50 b1, b2, b3 = self.recip_cell 

51 E_g = 0. 

52 for i in range(-N[0], N[0] + 1): 

53 for j in range(-N[1], N[1] + 1): 

54 for k in range(-N[2], N[2] + 1): 

55 g_v = 2 * pi * (i * b1 + j * b2 + k * b3) 

56 g2 = np.dot(g_v, g_v) 

57 if g2 > eps: # exclude g=0 

58 x = .25 * g2 / self.G**2 

59 E_g += np.exp(1.j * np.dot(r_v, g_v) - x) / x 

60 return pi / (self.Vcell * self.G**2) * E_g.real 

61 

62 def get_sum_real_ij(self, r_v, eps=1e-5): 

63 r"""The real space sum. 

64 

65 :: 

66 

67 ----- 

68 \ erfc( G [l-r| ) 

69 | -------------- 

70 / |l-r| 

71 ----- 

72 l not 0 

73 

74 Note: Add the l=0 term with erfc(r). 

75 """ 

76 N = self.Nl 

77 a1, a2, a3 = self.cell 

78 E_r = 0. 

79 for i in range(-N[0], N[0] + 1): 

80 for j in range(-N[1], N[1] + 1): 

81 for k in range(-N[2], N[2] + 1): 

82 l_v = i * a1 + j * a2 + k * a3 

83 if np.linalg.norm(l_v) > eps: # exclude l=0 

84 lr = np.linalg.norm(l_v - r_v) 

85 E_r += erfc(self.G * lr) / lr 

86 return E_r 

87 

88 def get_electrostatic_potential(self, r, r_B, q_B, excludefroml0=None): 

89 r"""... 

90 

91 Calculates the electrostatic potential at point r_i from point 

92 charges at {r_B} in a lattice using the Ewald summation. 

93 

94 Charge neutrality is obtained by adding the homogenous density 

95 q_hom/V:: 

96 

97 ---- ----' - 

98 \ \ q_j q_hom / 1 

99 phi(r_i) = | | --------------- + ----- |dr --------- 

100 / / |r_i - r_j + l| V | |r - r_i| 

101 ---- ---- / 

102 j l 

103 

104 r_B : matrix with the lattice basis (in cartesian coordinates). 

105 

106 q_B : probe charges (in units of e). 

107 

108 excludefroml0 : integer specifying if a point charge is not to 

109 be included in the central (l=0) unit cell. Used for Madelung 

110 constants. 

111 """ 

112 E0 = 0. 

113 if excludefroml0 is None: 

114 excludefroml0 = np.zeros(len(q_B), dtype=int) 

115 if excludefroml0 in range(len(q_B)): 

116 i = excludefroml0 

117 excludefroml0 = np.zeros(len(q_B), dtype=int) 

118 excludefroml0[i] = 1 

119 

120 assert sum(excludefroml0) <= 1 

121 

122 for i, q in enumerate(q_B): # potential from point charges 

123 rprime = r - r_B[i] 

124 absr = np.linalg.norm(rprime) 

125 E0 += q * self.get_sum_real_ij(rprime) 

126 E0 += q * self.get_sum_recip_ij(rprime) 

127 if excludefroml0[i]: # if sum over l not 0 

128 if absr < 1e-14: 

129 # lim r -> 0 erf(r G) / r = 2 * G / sqrt(pi) 

130 E0 -= q * 2. * self.G / np.sqrt(pi) 

131 else: 

132 E0 -= q * erf(absr * self.G) / absr 

133 else: # if sum over all l 

134 E0 += q * erfc(absr * self.G) / absr 

135 

136 # correct for compensating homogeneous background 

137 q_hom = -sum(q_B) 

138 E0 += q_hom * pi / (self.G**2 * self.Vcell) 

139 

140 return E0 

141 

142 

143def madelung(cell): 

144 return Ewald(cell).get_electrostatic_potential(np.zeros(3), np.zeros(3), 

145 [-1], 0)