Coverage for gpaw/pipekmezey/weightfunction.py: 97%

1r""" Class which builds a weight function for each atom 

2 in a molecule using simple atom centred gaussians. 

3 

4 This is the Hirshfeld scheme: 

5 A 

6 A p (r) 

7 w (r) = --------- 

8 mol 

9 \ n 

10 / p (r) 

11 -- 

12 n 

13 

14""" 

15 

16import numpy as np 

17from math import pi 

18from ase.units import Bohr 

19 

20# Cut-offs: [Ang] 

21Rc = {'Fe': 3.762, 

22 'O': 3.762, 

23 # 'H' : 3.762 

24 } 

25 

26# Gauss-width: [Ang] 

27mu = {'Fe': 1.0, 

28 'O': 0.6, 

29 # 'H' : 0.6 

30 } 

31 

32 

33class WeightFunc: 

34 

35 def __init__(self, gd, atoms, indexes, Rc=Rc, mu=mu): 

36 """ Given a grid-descriptor, atoms object and an index list 

37 construct a weight function defined by: 

38 

39 """ 

40 

41 self.gd = gd # Grid descriptor 

42 self.atoms = atoms 

43 self.indexes = indexes 

44 

45 # Construct Rc dict 

46 new = {} 

47 for a in self.atoms: 

48 if a.symbol in Rc: 

49 new[a.symbol] = Rc[a.symbol] 

50 else: 

51 new[a.symbol] = 3.762 

52 

53 self.Rc = new 

54 

55 # Construct mu dict 

56 new_mu = {} 

57 for a in self.atoms: 

58 if mu: 

59 if a.symbol in mu: 

60 new_mu[a.symbol] = mu[a.symbol] 

61 else: 

62 new_mu[a.symbol] = 0.85 

63 

64 # Larger atoms may need a bit more width? 

65 self.mu = new_mu 

66 

67 def truncated_gaussian(self, dis, mu, Rc): 

68 # Given mu and Rc construct Gaussian. 

69 # Gauss. is truncated at Rc 

70 

71 check = abs(dis) <= Rc / Bohr 

72 # Make gaussian 

73 gauss = 1.0 / (mu * np.sqrt(2.0 * pi)) * \ 

74 np.exp(- dis ** 2 / (2.0 * mu)) 

75 # Apply cut-off and send 

76 return (gauss * check) 

77 

78 def get_distance_vectors(self, pos): 

79 # Given atom position [Bohr], grab distances to all 

80 # grid-points (gpts) - employ MIC where appropriate. 

81 

82 # Scaled positions of gpts on some cpu, relative to all 

83 s_G = (np.indices(self.gd.n_c, float).T + 

84 self.gd.beg_c) / self.gd.N_c 

85 # print self.gd.n_c, self.gd.beg_c 

86 # Subtract scaled distance from atom to box boundary 

87 s_G -= np.linalg.solve(self.gd.cell_cv.T, pos) 

88 # MIC 

89 s_G -= self.gd.pbc_c * (2 * s_G).astype(int) 

90 # Apparently doing this twice works better... 

91 s_G -= self.gd.pbc_c * (2 * s_G).astype(int) 

92 # x,y,z distances 

93 xyz = np.dot(s_G, self.gd.cell_cv).T.copy() 

94 # 

95 return np.sqrt((xyz ** 2).sum(axis=0)) 

96 

97 def construct_total_density(self, atoms): 

98 # Add to empty grid 

99 empty = self.gd.zeros() 

100 

101 for atom in atoms: 

102 charge = atom.number 

103 symbol = atom.symbol 

104 

105 pos = atom.position / Bohr 

106 

107 dis = self.get_distance_vectors(pos) 

108 

109 empty += charge * self.truncated_gaussian( 

110 dis, self.mu[symbol], self.Rc[symbol]) 

111 # 

112 return empty 

113 

114 def construct_weight_function(self): 

115 # Grab atomic / molecular density 

116 dens_n = self.construct_total_density( 

117 self.atoms[self.indexes]) 

118 # Grab total density 

119 dens = self.construct_total_density(self.atoms) 

120 # Check zero elements 

121 check = dens == 0 

122 # Add constant to zeros ... 

123 dens += check * 1.0 

124 # make and send 

125 return dens_n / dens 

126 

127 

128class WignerSeitz: 

129 """ Construct weight function based on Wigner-Seitz 

130 

131 | 1 if (r-R) < (r-R') 

132 w(r) = | 

133 | 0 if (r-R) > (r-R') 

134 

135 for atom A at pos R 

136 

137 """ 

138 

139 def __init__(self, gd, atoms, index): 

140 # 

141 self.gd = gd 

142 self.atoms = atoms 

143 self.index = index 

144 

145 def get_distance_vectors(self, pos): 

146 # 

147 # Given atom position [Bohr], grab distances to all 

148 # grid-points (gpts) - employ MIC where appropriate. 

149 

150 # Scaled positions of gpts on some cpu, relative to all 

151 s_G = (np.indices(self.gd.n_c, float).T + 

152 self.gd.beg_c) / self.gd.N_c 

153 # 

154 # Subtract scaled distance from atom to box boundary 

155 s_G -= np.linalg.solve(self.gd.cell_cv.T, pos) 

156 # MIC 

157 s_G -= self.gd.pbc_c * (2 * s_G).astype(int) 

158 # Apparently doing this twice works better... 

159 s_G -= self.gd.pbc_c * (2 * s_G).astype(int) 

160 # x,y,z distances 

161 xyz = np.dot(s_G, self.gd.cell_cv).T.copy() 

162 # 

163 return np.sqrt((xyz ** 2).sum(axis=0)) 

164 

165 def construct_weight_function(self): 

166 # Grab distances of A to all gpts 

167 pos = self.atoms[self.index].position / Bohr 

168 dis_n = self.get_distance_vectors(pos) 

169 # 

170 empty = self.gd.zeros() 

171 norm = self.gd.zeros() 

172 # 

173 empty += 1 

174 norm += 1 

175 # 

176 for i, atom in enumerate(self.atoms): 

177 # 

178 if i == self.index: 

179 continue 

180 else: 

181 # 

182 pos = atom.position / Bohr 

183 dis_b = self.get_distance_vectors(pos) 

184 # 

185 check = dis_n <= dis_b # 

186 n_c = dis_n == dis_b # 

187 # 0 if longer, 1 if shorter, 1/no.(Ra=Ra') if same 

188 empty *= check 

189 norm += n_c * 1.0 

190 

191 return empty / norm