Coverage for gpaw/xc/tb09.py: 98%

1"""Tran-Blaha potential. 

2 

3From: 

4 

5 Accurate Band Gaps of Semiconductors and Insulators 

6 with a Semilocal Exchange-Correlation Potential 

7 

8 Fabien Tran and Peter Blaha 

9 

10 PRL 102, 226401 (2009) 

11 

12 DOI: 10.1103/PhysRevLett.102.226401 

13 

14""" 

15import numpy as np 

16from ase.utils import seterr 

17 

18from gpaw.xc.mgga import MGGA, weight_n 

19from gpaw.xc.libxc import LibXC 

20from gpaw.fd_operators import Laplace 

21 

22 

23class TB09Kernel: 

24 name = 'TB09' 

25 type = 'MGGA' 

26 alpha = -0.012 

27 beta = 1.023 

28 

29 def __init__(self, c=None): 

30 self.tb09 = LibXC('MGGA_X_TB09', provides_laplacian=True).xc.tb09 

31 self.ldac = LibXC('LDA_C_PW') 

32 

33 self.fixedc = c is not None # calculate c or use fixed value 

34 self.c = c # amount of "exact exchange" 

35 self.n = 0 # Lebedev quadrature point number (0-49) 

36 self.sign = 1.0 # sign of PAW correction: +1 for AE and -1 for PS 

37 self.I = None # integral from Eq. (3) 

38 

39 def calculate(self, e_g, n_sg, dedn_sg, sigma_xg, 

40 dedsigma_xg, tau_sg, dedtau_sg): 

41 ns = len(n_sg) 

42 n_sg[n_sg < 1e-6] = 1e-6 

43 

44 if n_sg.ndim == 4: 

45 if not self.fixedc: 

46 if self.c is None: 

47 # We don't have the integral yet - just use 1.0: 

48 self.c = 1.0 

49 else: 

50 self.I = self.world.sum_scalar(self.I) 

51 self.c = (self.alpha + self.beta * 

52 (self.I / self.gd.volume)**0.5) 

53 

54 # Start calculation of c for use in the next SCF step: 

55 if ns == 1: 

56 gradn_g = sigma_xg[0]**0.5 

57 else: 

58 gradn_g = (sigma_xg[0] + 

59 2 * sigma_xg[1] + 

60 sigma_xg[2])**0.5 

61 self.I = self.gd.integrate(gradn_g / n_sg.sum(0)) 

62 # The domain is not distributed like the PAW corrections: 

63 self.I /= self.world.size 

64 

65 lapl_sg = self.gd.empty(ns) 

66 for n_g, lapl_g in zip(n_sg, lapl_sg): 

67 self.lapl.apply(n_g, lapl_g) 

68 

69 else: 

70 rgd = self.rgd 

71 lapl_sg = [] 

72 for n_Lg in self.n_sLg: 

73 lapl_g = rgd.laplace(np.dot(self.Y_L, n_Lg)) 

74 l = 0 

75 L1 = 0 

76 while L1 < len(self.Y_L): 

77 L2 = L1 + 2 * l + 1 

78 n_g = np.dot(self.Y_L[L1:L2], n_Lg[L1:L2]) 

79 with seterr(divide='ignore', invalid='ignore'): 

80 lapl_g -= l * (l + 1) * n_g / rgd.r_g**2 

81 lapl_g[0] = 0.0 

82 L1 = L2 

83 l += 1 

84 lapl_sg.append(lapl_g) 

85 

86 if not self.fixedc: 

87 # PAW corrections to integral: 

88 w = self.sign * weight_n[self.n] 

89 if ns == 1: 

90 gradn_g = sigma_xg[0]**0.5 

91 else: 

92 gradn_g = (sigma_xg[0] + 

93 2 * sigma_xg[1] + 

94 sigma_xg[2])**0.5 

95 self.I += w * rgd.integrate(gradn_g / n_sg.sum(0)) 

96 

97 self.n += 1 

98 if self.n == len(weight_n): 

99 self.n = 0 

100 self.sign = -self.sign 

101 

102 # dedn_sg[:] = 0.0 

103 sigma_xg[sigma_xg < 1e-10] = 1e-10 

104 tau_sg[tau_sg < 1e-10] = 1e-10 

105 

106 for n_g, sigma_g, lapl_g, tau_g, v_g in zip(n_sg, sigma_xg[::2], 

107 lapl_sg, tau_sg, dedn_sg): 

108 self.tb09(self.c, n_g.ravel(), sigma_g, lapl_g, tau_g, v_g, 

109 dedsigma_xg) 

110 

111 self.ldac.calculate(e_g, n_sg, dedn_sg) 

112 e_g[:] = 0.0 

113 

114 dedsigma_xg[:] = 0.0 

115 dedtau_sg[:] = 0.0 

116 

117 

118class TB09(MGGA): 

119 def __init__(self, c=None, stencil=2): 

120 MGGA.__init__(self, TB09Kernel(c), stencil=stencil) 

121 

122 def get_setup_name(self): 

123 return 'LDA' 

124 

125 def initialize(self, dens, ham, wfs): 

126 MGGA.initialize(self, dens, ham, wfs) 

127 self.kernel.world = wfs.world 

128 self.kernel.gd = dens.finegd 

129 self.kernel.lapl = Laplace(dens.finegd) 

130 

131 def create_mgga_radial_calculator(self): 

132 rcalc = MGGA.create_mgga_radial_calculator(self) 

133 

134 def f(rgd, n_sLg, Y_L, dndr_sLg, rnablaY_Lv, n): 

135 self.kernel.n_sLg = n_sLg 

136 self.kernel.Y_L = Y_L 

137 self.kernel.rgd = rgd 

138 return rcalc(rgd, n_sLg, Y_L, dndr_sLg, rnablaY_Lv, n) 

139 

140 return f 

141 

142 def apply_orbital_dependent_hamiltonian(self, kpt, psit_xG, 

143 Htpsit_xG, dH_asp): 

144 pass 

145 

146 @property 

147 def c(self): 

148 """Amount of "exact exchange".""" 

149 return self.kernel.c