Coverage for gpaw/xc/bee.py: 96%

1import numpy as np 

2from ase.units import Hartree 

3 

4import gpaw.cgpaw as cgpaw 

5from gpaw.xc import XC 

6from gpaw.xc.kernel import XCKernel 

7from gpaw.xc.libxc import LibXC 

8from gpaw.xc.vdw import VDWFunctional 

9from gpaw import debug 

10 

11 

12class BEE2(XCKernel): 

13 """GGA exchange expanded in Legendre polynomials.""" 

14 def __init__(self, parameters=None): 

15 """BEE2. 

16 

17 parameters: array 

18 [transformation,0.0,[orders],[coefs]]. 

19 

20 """ 

21 

22 if parameters is None: 

23 parameters = ( 

24 [4.0, 0.0] + list(range(30)) + 

25 [1.516501714304992365356, 0.441353209874497942611, 

26 -0.091821352411060291887, -0.023527543314744041314, 

27 0.034188284548603550816, 0.002411870075717384172, 

28 -0.014163813515916020766, 0.000697589558149178113, 

29 0.009859205136982565273, -0.006737855050935187551, 

30 -0.001573330824338589097, 0.005036146253345903309, 

31 -0.002569472452841069059, -0.000987495397608761146, 

32 0.002033722894696920677, -0.000801871884834044583, 

33 -0.000668807872347525591, 0.001030936331268264214, 

34 -0.000367383865990214423, -0.000421363539352619543, 

35 0.000576160799160517858, -0.000083465037349510408, 

36 -0.000445844758523195788, 0.000460129009232047457, 

37 -0.000005231775398304339, -0.000423957047149510404, 

38 0.000375019067938866537, 0.000021149381251344578, 

39 -0.000190491156503997170, 0.000073843624209823442]) 

40 else: 

41 assert len(parameters) > 2 

42 assert np.mod(len(parameters), 2) == 0 

43 assert parameters[1] == 0.0 

44 

45 parameters = np.array(parameters, dtype=float).ravel() 

46 self.xc = cgpaw.XCFunctional(17, parameters) 

47 self.type = 'GGA' 

48 self.name = 'BEE2' 

49 

50 

51class BEEVDWKernel(XCKernel): 

52 """Kernel for BEEVDW functionals.""" 

53 def __init__(self, bee, xcoefs, ldac, ggac): 

54 """BEEVDW kernel. 

55 

56 parameters: 

57 

58 bee : str 

59 choose BEE1 or BEE2 exchange basis expansion. 

60 xcoefs : array 

61 coefficients for exchange. 

62 ldac : float 

63 coefficient for LDA correlation. 

64 pbec : float 

65 coefficient for PBE correlation. 

66 

67 """ 

68 

69 if bee == 'BEE2': 

70 self.BEE = BEE2(xcoefs) 

71 self.GGAc = LibXC('GGA_C_PBE') 

72 self.xtype = 'GGA' 

73 self.type = 'GGA' 

74 elif bee == 'BEE3': 

75 self.BEE = LibXC('MGGA_X_MBEEFVDW') 

76 self.GGAc = LibXC('GGA_C_PBE_SOL') 

77 self.xtype = 'MGGA' 

78 self.type = 'MGGA' 

79 else: 

80 raise ValueError('Unknown BEE exchange: %s', bee) 

81 

82 self.LDAc = LibXC('LDA_C_PW') 

83 self.ldac = ldac 

84 self.ggac = ggac 

85 if bee in ['BEE1', 'BEE2']: 

86 self.ggac -= 1.0 

87 self.name = 'BEEVDW' 

88 

89 def calculate(self, e_g, n_sg, dedn_sg, 

90 sigma_xg=None, dedsigma_xg=None, 

91 tau_sg=None, dedtau_sg=None): 

92 if debug: 

93 self.check_arguments(e_g, n_sg, dedn_sg, sigma_xg, dedsigma_xg, 

94 tau_sg, dedtau_sg) 

95 

96 if self.xtype == 'GGA': 

97 self.BEE.calculate(e_g, n_sg, dedn_sg, sigma_xg, dedsigma_xg) 

98 elif self.xtype == 'MGGA': 

99 self.BEE.calculate(e_g, n_sg, dedn_sg, sigma_xg, dedsigma_xg, 

100 tau_sg, dedtau_sg) 

101 else: 

102 raise ValueError('Unexpected value of xtype:', self.xtype) 

103 

104 e0_g = np.empty_like(e_g) 

105 dedn0_sg = np.empty_like(dedn_sg) 

106 dedsigma0_xg = np.empty_like(dedsigma_xg) 

107 for coef, kernel in [(self.ldac, self.LDAc), 

108 (self.ggac, self.GGAc)]: 

109 dedn0_sg[:] = 0.0 

110 kernel.calculate(e0_g, n_sg, dedn0_sg, sigma_xg, dedsigma0_xg) 

111 e_g += coef * e0_g 

112 dedn_sg += coef * dedn0_sg 

113 if kernel.type == 'GGA': 

114 dedsigma_xg += coef * dedsigma0_xg 

115 

116 

117class BEEFEnsemble: 

118 """BEEF ensemble error estimation.""" 

119 def __init__(self, calc): 

120 """BEEF ensemble 

121 

122 parameters: 

123 

124 calc: object 

125 Calculator holding a selfconsistent BEEF type electron density. 

126 May be BEEF-vdW or mBEEF. 

127 """ 

128 

129 self.calc = calc 

130 

131 self.e_dft = None 

132 self.e0 = None 

133 

134 # determine functional and read parameters 

135 self.xc = self.calc.get_xc_functional() 

136 if self.xc == 'BEEF-vdW': 

137 self.bee_type = 1 

138 elif self.xc == 'mBEEF': 

139 self.bee_type = 2 

140 self.max_order = 8 

141 self.trans = [6.5124, -1.0] 

142 self.calc.converge_wave_functions() 

143 self.calc.log('wave functions converged') 

144 elif self.xc == 'mBEEF-vdW': 

145 self.bee_type = 3 

146 self.max_order = 5 

147 self.trans = [6.5124, -1.0] 

148 self.calc.converge_wave_functions() 

149 self.calc.log('wave functions converged') 

150 else: 

151 raise NotImplementedError('xc = %s not implemented' % self.xc) 

152 

153 def create_xc_contributions(self, type): 

154 """General function for creating exchange or correlation energies""" 

155 assert type in ['exch', 'corr'] 

156 err = 'bee_type %i not implemented' % self.bee_type 

157 

158 if type == 'exch': 

159 if self.bee_type == 1: 

160 out = self.beefvdw_energy_contribs_x() 

161 elif self.bee_type in [2, 3]: 

162 out = self.mbeef_exchange_energy_contribs() 

163 else: 

164 raise NotImplementedError(err) 

165 else: 

166 if self.bee_type == 1: 

167 out = self.beefvdw_energy_contribs_c() 

168 elif self.bee_type == 2: 

169 out = np.array([]) 

170 elif self.bee_type == 3: 

171 out = self.mbeefvdw_energy_contribs_c() 

172 else: 

173 raise NotImplementedError(err) 

174 return out 

175 

176 def get_non_xc_total_energies(self): 

177 """Compile non-XC total energy contributions""" 

178 if self.e_dft is None: 

179 self.e_dft = self.calc.get_potential_energy() 

180 if self.e0 is None: 

181 from gpaw.xc.kernel import XCNull 

182 xc_null = XC(XCNull()) 

183 self.e0 = self.e_dft + self.calc.get_xc_difference(xc_null) 

184 assert isinstance(self.e_dft, float) 

185 assert isinstance(self.e0, float) 

186 

187 def mbeef_exchange_energy_contribs(self): 

188 """Legendre polynomial exchange contributions to mBEEF Etot""" 

189 self.get_non_xc_total_energies() 

190 e_x = np.zeros((self.max_order, self.max_order)) 

191 for p1 in range(self.max_order): # alpha 

192 for p2 in range(self.max_order): # s2 

193 pars_i = np.array([1, self.trans[0], p2, 1.0]) 

194 pars_j = np.array([1, self.trans[1], p1, 1.0]) 

195 pars = np.hstack((pars_i, pars_j)) 

196 x = XC('2D-MGGA', pars) 

197 e_x[p1, p2] = (self.e_dft + 

198 self.calc.get_xc_difference(x) - self.e0) 

199 del x 

200 return e_x 

201 

202 def beefvdw_energy_contribs_x(self): 

203 """Legendre polynomial exchange contributions to BEEF-vdW Etot""" 

204 self.get_non_xc_total_energies() 

205 e_pbe = (self.e_dft + self.calc.get_xc_difference('GGA_C_PBE') - 

206 self.e0) 

207 

208 exch = np.zeros(30) 

209 for p in range(30): 

210 pars = [4, 0, p, 1.0] 

211 bee = XC('BEE2', pars) 

212 exch[p] = (self.e_dft + self.calc.get_xc_difference(bee) - 

213 self.e0 - e_pbe) 

214 del bee 

215 return exch 

216 

217 def beefvdw_energy_contribs_c(self): 

218 """LDA and PBE correlation contributions to BEEF-vdW Etot""" 

219 self.get_non_xc_total_energies() 

220 e_lda = self.e_dft + self.calc.get_xc_difference('LDA_C_PW') - self.e0 

221 e_pbe = self.e_dft + self.calc.get_xc_difference('GGA_C_PBE') - self.e0 

222 corr = np.array([e_lda, e_pbe]) 

223 return corr 

224 

225 def mbeefvdw_energy_contribs_c(self): 

226 """LDA, PBEsol, and nl2 correlation contributions to mBEEF-vdW Etot""" 

227 self.get_non_xc_total_energies() 

228 e_lda = self.e_dft + self.calc.get_xc_difference('LDA_C_PW') - self.e0 

229 e_sol = (self.e_dft + self.calc.get_xc_difference('GGA_C_PBE_SOL') - 

230 self.e0) 

231 vdwdf2 = VDWFunctional('vdW-DF2') 

232 self.calc.get_xc_difference(vdwdf2) 

233 e_nl2 = vdwdf2.get_Ecnl() * Hartree 

234 corr = np.array([e_lda, e_sol, e_nl2]) 

235 return corr