Coverage for gpaw/response/temp.py: 99%

1import numpy as np 

2from functools import cached_property 

3from gpaw.response.pw_parallelization import Blocks1D 

4from gpaw.response.gamma_int import GammaIntegral 

5 

6 

7class DielectricFunctionCalculator: 

8 def __init__(self, chi0, coulomb, xckernel, mode): 

9 self.coulomb = coulomb 

10 self.qpd = chi0.qpd 

11 self.mode = mode 

12 self.optical_limit = chi0.optical_limit 

13 self.chi0 = chi0 

14 self.xckernel = xckernel 

15 self.wblocks = Blocks1D(chi0.body.blockdist.blockcomm, len(chi0.wd)) 

16 # Generate fine grid in vicinity of gamma 

17 if chi0.optical_limit and self.wblocks.nlocal: 

18 self.gamma_int = GammaIntegral(self.coulomb, self.qpd) 

19 else: 

20 self.gamma_int = None 

21 

22 @cached_property 

23 def sqrtV_G(self): 

24 return self.coulomb.sqrtV(qpd=self.qpd, q_v=None) 

25 

26 @cached_property 

27 def I_GG(self): 

28 return np.eye(self.qpd.ngmax) 

29 

30 @cached_property 

31 def fxc_GG(self): 

32 if self.mode == 'GW': 

33 return self.I_GG 

34 else: 

35 return self.xckernel.calculate(self.qpd) 

36 

37 @cached_property 

38 def chi0_wGG(self): 

39 return self.chi0.body.copy_array_with_distribution('wGG') 

40 

41 def get_epsinv_wGG(self, only_correlation=True): 

42 """ 

43 Calculates inverse dielectric matrix for all frequencies. 

44 """ 

45 epsinv_wGG = [] 

46 for w in range(self.wblocks.nlocal): 

47 epsinv_GG = self.single_frequency_epsinv_GG(w) 

48 if only_correlation: 

49 epsinv_GG -= self.I_GG 

50 epsinv_wGG.append(epsinv_GG) 

51 return np.asarray(epsinv_wGG) 

52 

53 def single_frequency_epsinv_GG(self, w): 

54 """ 

55 Calculates inverse dielectric matrix for single frequency 

56 """ 

57 _dfc = _DielectricFunctionCalculator(self.sqrtV_G, 

58 self.chi0_wGG[w], 

59 self.mode, 

60 self.fxc_GG) 

61 if self.optical_limit: 

62 W = self.wblocks.a + w 

63 _dfc = _GammaDielectricFunctionCalculator( 

64 _dfc, self.gamma_int, 

65 self.chi0.chi0_Wvv[W], self.chi0.chi0_WxvG[W]) 

66 return _dfc.get_epsinv_GG() 

67 

68 

69class _DielectricFunctionCalculator: 

70 def __init__(self, sqrtV_G, chi0_GG, mode, fxc_GG=None): 

71 self.sqrtV_G = sqrtV_G 

72 self.chiVV_GG = chi0_GG * sqrtV_G * sqrtV_G[:, np.newaxis] 

73 

74 self.I_GG = np.eye(len(sqrtV_G)) 

75 

76 self.fxc_GG = fxc_GG 

77 self.chi0_GG = chi0_GG 

78 self.mode = mode 

79 

80 def new_with(self, *, sqrtV_G, chi0_GG): 

81 return _DielectricFunctionCalculator( 

82 sqrtV_G, chi0_GG, self.mode, fxc_GG=self.fxc_GG) 

83 

84 def _chiVVfxc_GG(self): 

85 assert self.mode != 'GW' 

86 assert self.fxc_GG is not None 

87 return self.chiVV_GG @ self.fxc_GG 

88 

89 def eps_GG_gwp(self): 

90 gwp_inv_GG = np.linalg.inv(self.I_GG - self._chiVVfxc_GG() + 

91 self.chiVV_GG) 

92 return self.I_GG - gwp_inv_GG @ self.chiVV_GG 

93 

94 def eps_GG_gws(self): 

95 # Note how the signs are different wrt. gwp. 

96 # Nobody knows why. 

97 gws_inv_GG = np.linalg.inv(self.I_GG + self._chiVVfxc_GG() - 

98 self.chiVV_GG) 

99 return gws_inv_GG @ (self.I_GG - self.chiVV_GG) 

100 

101 def eps_GG_plain(self): 

102 return self.I_GG - self.chiVV_GG 

103 

104 def eps_GG_w_fxc(self): 

105 return self.I_GG - self._chiVVfxc_GG() 

106 

107 def get_eps_GG(self): 

108 mode = self.mode 

109 if mode == 'GWP': 

110 return self.eps_GG_gwp() 

111 elif mode == 'GWS': 

112 return self.eps_GG_gws() 

113 elif mode == 'GW': 

114 return self.eps_GG_plain() 

115 elif mode == 'GWG': 

116 return self.eps_GG_w_fxc() 

117 raise ValueError(f'Unknown mode: {mode}') 

118 

119 def get_epsinv_GG(self): 

120 eps_GG = self.get_eps_GG() 

121 return np.linalg.inv(eps_GG) 

122 

123 

124class _GammaDielectricFunctionCalculator: 

125 

126 def __init__(self, _dfc, gamma_int, chi0_vv, chi0_xvG): 

127 self._dfc = _dfc 

128 self.gamma_int = gamma_int 

129 

130 self.chi0_vv = chi0_vv 

131 self.chi0_xvG = chi0_xvG 

132 

133 @property 

134 def chi0_GG(self): 

135 return self._dfc.chi0_GG 

136 

137 def get_epsinv_GG(self): 

138 # Get average epsinv over small region around Gamma 

139 epsinv_GG = np.zeros(self.chi0_GG.shape, complex) 

140 for qweight, sqrtV_G, chi0_mapping in self.gamma_int: 

141 chi0p_GG = chi0_mapping(self.chi0_GG, self.chi0_vv, self.chi0_xvG) 

142 _dfc = self._dfc.new_with(sqrtV_G=sqrtV_G, chi0_GG=chi0p_GG) 

143 epsinv_GG += qweight * _dfc.get_epsinv_GG() 

144 return epsinv_GG