Coverage for gpaw/xc/gllb/c_gllbscr.py: 99%

1import numpy as np 

2from math import sqrt, pi 

3 

4from gpaw.xc.pawcorrection import rnablaY_nLv 

5from gpaw.sphere.lebedev import weight_n 

6from gpaw.xc import XC 

7from gpaw.xc.gllb.contribution import Contribution 

8 

9 

10class C_GLLBScr(Contribution): 

11 def __init__(self, weight, functional, damp=1e-10): 

12 Contribution.__init__(self, weight) 

13 self.xc = XC(functional) 

14 self.damp = damp 

15 

16 def set_damp(self, damp): 

17 self.damp = damp 

18 

19 def get_name(self): 

20 return 'SCREENING' 

21 

22 def get_desc(self): 

23 desc = self.xc.get_description() 

24 return self.xc.name if desc is None else desc 

25 

26 def initialize(self, density, hamiltonian, wfs): 

27 Contribution.initialize(self, density, hamiltonian, wfs) 

28 # Always 1 spin, no matter what calculation nspins is 

29 self.vt_sg = self.finegd.empty(1) 

30 self.e_g = self.finegd.empty() 

31 

32 def initialize_1d(self, ae): 

33 Contribution.initialize_1d(self, ae) 

34 self.v_g = np.zeros(self.ae.N) 

35 self.e_g = np.zeros(self.ae.N) 

36 

37 # Calcualte the GLLB potential and energy 1d 

38 def add_xc_potential_and_energy_1d(self, v_g): 

39 self.v_g[:] = 0.0 

40 self.e_g[:] = 0.0 

41 self.xc.calculate_spherical(self.ae.rgd, self.ae.n.reshape((1, -1)), 

42 self.v_g.reshape((1, -1)), self.e_g) 

43 v_g += 2 * self.weight * self.e_g / (self.ae.n + self.damp) 

44 Exc = self.weight * np.sum(self.e_g * self.ae.rgd.dv_g) 

45 return Exc 

46 

47 def calculate(self, e_g, n_sg, v_sg): 

48 # Calculate spin-paired exchange screening 

49 # as Eqs. (21-22) of https://doi.org/10.1103/PhysRevB.82.115106 

50 # and spin-polarized exchange screening 

51 # as two spin-paired calculations n=2*n_s 

52 mult = self.nspins # mult = 1 for spin-paired and 2 for spin-polarized 

53 for n_g, v_g in zip(n_sg, v_sg): 

54 self.e_g[:] = 0.0 

55 self.vt_sg[:] = 0.0 # Note: this array has nspins=1 always 

56 self.xc.calculate(self.finegd, mult * n_g[np.newaxis], 

57 self.vt_sg, self.e_g) 

58 self.e_g[:] = np.where(n_g < self.damp, 0, self.e_g) 

59 v_g += self.weight * 2 * self.e_g / (mult * n_g + self.damp) 

60 e_g += self.weight * self.e_g / mult 

61 

62 def calculate_energy_and_derivatives(self, setup, D_sp, H_sp, a, 

63 addcoredensity=True): 

64 # Get the XC-correction instance 

65 c = setup.xc_correction 

66 nspins = self.nspins 

67 

68 E = 0 

69 for D_p, dEdD_p in zip(D_sp, H_sp): 

70 D_Lq = np.dot(c.B_pqL.T, nspins * D_p) 

71 n_Lg = np.dot(D_Lq, c.n_qg) 

72 if addcoredensity: 

73 n_Lg[0] += c.nc_g * sqrt(4 * pi) 

74 nt_Lg = np.dot(D_Lq, c.nt_qg) 

75 if addcoredensity: 

76 nt_Lg[0] += c.nct_g * sqrt(4 * pi) 

77 dndr_Lg = np.zeros((c.Lmax, c.ng)) 

78 dntdr_Lg = np.zeros((c.Lmax, c.ng)) 

79 for L in range(c.Lmax): 

80 c.rgd.derivative(n_Lg[L], dndr_Lg[L]) 

81 c.rgd.derivative(nt_Lg[L], dntdr_Lg[L]) 

82 vt_g = np.zeros(c.ng) 

83 v_g = np.zeros(c.ng) 

84 e_g = np.zeros(c.ng) 

85 deda2_g = np.zeros(c.ng) 

86 for y, (w, Y_L) in enumerate(zip(weight_n, c.Y_nL)): 

87 # Cut gradient releated coefficient to match the setup's Lmax 

88 A_Li = rnablaY_nLv[y, :c.Lmax] 

89 

90 # Expand pseudo density 

91 nt_g = np.dot(Y_L, nt_Lg) 

92 

93 # Expand pseudo density gradient 

94 a1x_g = np.dot(A_Li[:, 0], nt_Lg) 

95 a1y_g = np.dot(A_Li[:, 1], nt_Lg) 

96 a1z_g = np.dot(A_Li[:, 2], nt_Lg) 

97 a2_g = a1x_g**2 + a1y_g**2 + a1z_g**2 

98 a2_g[1:] /= c.rgd.r_g[1:]**2 

99 a2_g[0] = a2_g[1] 

100 a1_g = np.dot(Y_L, dntdr_Lg) 

101 a2_g += a1_g**2 

102 

103 vt_g[:] = 0.0 

104 e_g[:] = 0.0 

105 # Calculate pseudo GGA energy density (potential is discarded) 

106 self.xc.kernel.calculate(e_g, nt_g.reshape((1, -1)), 

107 vt_g.reshape((1, -1)), 

108 a2_g.reshape((1, -1)), 

109 deda2_g.reshape((1, -1))) 

110 

111 # Calculate pseudo GLLB-potential from GGA-energy density 

112 vt_g[:] = 2 * e_g / (nt_g + self.damp) 

113 

114 dEdD_p -= self.weight * w * np.dot(np.dot(c.B_pqL, Y_L), 

115 np.dot(c.nt_qg, 

116 vt_g * c.rgd.dv_g)) 

117 

118 E -= w * np.dot(e_g, c.rgd.dv_g) / nspins 

119 

120 # Expand density 

121 n_g = np.dot(Y_L, n_Lg) 

122 

123 # Expand density gradient 

124 a1x_g = np.dot(A_Li[:, 0], n_Lg) 

125 a1y_g = np.dot(A_Li[:, 1], n_Lg) 

126 a1z_g = np.dot(A_Li[:, 2], n_Lg) 

127 a2_g = a1x_g**2 + a1y_g**2 + a1z_g**2 

128 a2_g[1:] /= c.rgd.r_g[1:]**2 

129 a2_g[0] = a2_g[1] 

130 a1_g = np.dot(Y_L, dndr_Lg) 

131 a2_g += a1_g**2 

132 

133 v_g[:] = 0.0 

134 e_g[:] = 0.0 

135 # Calculate GGA energy density (potential is discarded) 

136 self.xc.kernel.calculate(e_g, n_g.reshape((1, -1)), 

137 v_g.reshape((1, -1)), 

138 a2_g.reshape((1, -1)), 

139 deda2_g.reshape((1, -1))) 

140 

141 # Calculate GLLB-potential from GGA-energy density 

142 v_g[:] = 2 * e_g / (n_g + self.damp) 

143 

144 dEdD_p += self.weight * w * np.dot(np.dot(c.B_pqL, Y_L), 

145 np.dot(c.n_qg, 

146 v_g * c.rgd.dv_g)) 

147 E += w * np.dot(e_g, c.rgd.dv_g) / nspins 

148 

149 return E * self.weight 

150 

151 def add_smooth_xc_potential_and_energy_1d(self, vt_g): 

152 self.v_g[:] = 0.0 

153 self.e_g[:] = 0.0 

154 self.xc.calculate_spherical(self.ae.rgd, self.ae.nt.reshape((1, -1)), 

155 self.v_g.reshape((1, -1)), self.e_g) 

156 vt_g += 2 * self.weight * self.e_g / (self.ae.nt + self.damp) 

157 return self.weight * np.sum(self.e_g * self.ae.rgd.dv_g)