Coverage for gpaw/response/density_kernels.py: 93%

1"""XC density kernels for response function calculations.""" 

2from abc import ABC, abstractmethod 

3from dataclasses import dataclass 

4 

5import numpy as np 

6 

7from gpaw.response import ResponseGroundStateAdapter, ResponseContext 

8from gpaw.response.qpd import SingleQPWDescriptor 

9from gpaw.response.localft import LocalFTCalculator 

10from gpaw.response.fxc_kernels import AdiabaticFXCCalculator 

11 

12 

13@dataclass 

14class DensityXCKernel(ABC): 

15 gs: ResponseGroundStateAdapter 

16 context: ResponseContext 

17 functional: str 

18 

19 def __post_init__(self): 

20 assert self.gs.nspins == 1 

21 

22 @staticmethod 

23 def from_functional(gs, context, functional, **kwargs): 

24 """Factory function creating DensityXCKernels. 

25 

26 Choose between ALDA, Bootstrap and LRalpha (long-range kernel), where 

27 alpha is a user specified parameter (for example functional='LR0.25'). 

28 """ 

29 if functional[0] == 'A': 

30 return AdiabaticDensityKernel(gs, context, functional, kwargs) 

31 elif functional[:2] == 'LR': 

32 return LRDensityKernel(gs, context, functional) 

33 elif functional == 'Bootstrap': 

34 return BootstrapDensityKernel(gs, context, functional) 

35 raise ValueError( 

36 'Invalid functional for the density-density xc kernel:' 

37 f'{functional}') 

38 

39 def __call__(self, qpd: SingleQPWDescriptor, chi0_wGG=None): 

40 self.context.print(f'Calculating {self}') 

41 return self.calculate(qpd=qpd, chi0_wGG=chi0_wGG) 

42 

43 @abstractmethod 

44 def __repr__(self): 

45 """String representation of the density xc kernel.""" 

46 

47 @abstractmethod 

48 def calculate(self, *args, **kwargs): 

49 """Calculate the xc kernel. 

50 

51 Since the exchange-correlation kernel is going to be rescaled according 

52 to the bare Coulomb interaction, 

53 ˷ 

54 K_xc(q) = v^(-1/2)(q) K_xc(q) v^(-1/2)(q) 

55 

56 and that the Coulomb interaction in the optical q→0 limit leaves the 

57 long-range q-dependence out, see gpaw.response.coulomb_kernels, 

58 

59 v(q) = 4π/|G+q| -> 4π for G==0 and 4π/|G| otherwise if q==0, 

60 

61 we need to account to account for it here. That is, 

62 

63 For q→0, the head and wings of the returned K_xc(q) needs to be 

64 rescaled according to K_xc(q) -> q K_xc(q) q 

65 """ 

66 

67 

68@dataclass 

69class AdiabaticDensityKernel(DensityXCKernel): 

70 rshe_kwargs: dict 

71 

72 def __post_init__(self): 

73 super().__post_init__() 

74 localft_calc = LocalFTCalculator.from_rshe_parameters( 

75 self.gs, self.context, **self.rshe_kwargs) 

76 self.fxc_calc = AdiabaticFXCCalculator(localft_calc) 

77 

78 def __repr__(self): 

79 return f'{self.functional} kernel' 

80 

81 def calculate(self, qpd: SingleQPWDescriptor, **ignored): 

82 fxc_kernel = self.fxc_calc(self.functional, '00', qpd) 

83 Kxc_GG = fxc_kernel.get_Kxc_GG() 

84 if qpd.optical_limit: 

85 Kxc_GG[0, :] = 0.0 

86 Kxc_GG[:, 0] = 0.0 

87 return Kxc_GG 

88 

89 

90@dataclass 

91class LRDensityKernel(DensityXCKernel): 

92 

93 def __post_init__(self): 

94 super().__post_init__() 

95 self.alpha = float(self.functional[2:]) 

96 

97 def __repr__(self): 

98 return f'LR kernel with alpha = {self.alpha}' 

99 

100 def calculate(self, qpd: SingleQPWDescriptor, **ignored): 

101 Kxc_sGG = calculate_lr_kernel(qpd, alpha=self.alpha) 

102 return Kxc_sGG[0] 

103 

104 

105@dataclass 

106class BootstrapDensityKernel(DensityXCKernel): 

107 

108 def __repr__(self): 

109 return 'Bootstrap kernel' 

110 

111 def calculate(self, qpd, *, chi0_wGG): 

112 Kxc_sGG = get_bootstrap_kernel(qpd, chi0_wGG, self.context) 

113 return Kxc_sGG[0] 

114 

115 

116def calculate_lr_kernel(qpd, alpha=0.2): 

117 """Long range kernel: fxc = \alpha / |q+G|^2""" 

118 

119 assert qpd.optical_limit 

120 

121 f_G = np.zeros(len(qpd.G2_qG[0])) 

122 f_G[0] = -alpha 

123 f_G[1:] = -alpha / qpd.G2_qG[0][1:] 

124 

125 return np.array([np.diag(f_G)]) 

126 

127 

128def get_bootstrap_kernel(qpd, chi0_wGG, context): 

129 """ Bootstrap kernel (see below) """ 

130 

131 if context.comm.rank == 0: 

132 chi0_GG = chi0_wGG[0] 

133 if context.comm.size > 1: 

134 # If size == 1, chi0_GG is not contiguous, and broadcast() 

135 # will fail in debug mode. So we skip it until someone 

136 # takes a closer look. 

137 context.comm.broadcast(chi0_GG, 0) 

138 else: 

139 nG = qpd.ngmax 

140 chi0_GG = np.zeros((nG, nG), complex) 

141 context.comm.broadcast(chi0_GG, 0) 

142 

143 return calculate_bootstrap_kernel(qpd, chi0_GG, context) 

144 

145 

146def calculate_bootstrap_kernel(qpd, chi0_GG, context): 

147 """Bootstrap kernel PRL 107, 186401""" 

148 p = context.print 

149 

150 if qpd.optical_limit: 

151 v_G = np.zeros(len(qpd.G2_qG[0])) 

152 v_G[0] = 4 * np.pi 

153 v_G[1:] = 4 * np.pi / qpd.G2_qG[0][1:] 

154 else: 

155 v_G = 4 * np.pi / qpd.G2_qG[0] 

156 

157 nG = len(v_G) 

158 K_GG = np.diag(v_G) 

159 

160 Kxc_GG = np.zeros((nG, nG), dtype=complex) 

161 dminv_GG = np.zeros((nG, nG), dtype=complex) 

162 

163 for iscf in range(120): 

164 dminvold_GG = dminv_GG.copy() 

165 Kxc_GG = K_GG + Kxc_GG 

166 

167 chi_GG = np.dot(np.linalg.inv(np.eye(nG, nG) 

168 - np.dot(chi0_GG, Kxc_GG)), chi0_GG) 

169 dminv_GG = np.eye(nG, nG) + np.dot(K_GG, chi_GG) 

170 

171 alpha = dminv_GG[0, 0] / (K_GG[0, 0] * chi0_GG[0, 0]) 

172 Kxc_GG = alpha * K_GG 

173 p(iscf, 'alpha =', alpha, flush=False) 

174 error = np.abs(dminvold_GG - dminv_GG).sum() 

175 if np.sum(error) < 0.1: 

176 p('Self consistent fxc finished in %d iterations !' % iscf) 

177 break 

178 if iscf > 100: 

179 p('Too many fxc scf steps !') 

180 

181 return np.array([Kxc_GG])