Coverage for gpaw/xc/gllb/coefficients.py: 89%

1from __future__ import annotations 

2 

3import numpy as np 

4from ase.units import Ha 

5from scipy.special import erfcx 

6 

7K_G = 8 * np.sqrt(2) / (3 * np.pi**2) # 0.382106112167171 

8 

9 

10class Coefficients: 

11 """Coefficient calculator for GLLB functionals. 

12 

13 This class implements the calculation of sqrt(E) coefficients as given by 

14 Eq. (16) of https://doi.org/10.1103/PhysRevB.82.115106. 

15 

16 Parameters 

17 ---------- 

18 eps (in eV) 

19 This parameter cuts sqrt(E) to zero for E < eps. 

20 The cut is supposed to help convergence with degenerate systems. 

21 This parameter should be small. 

22 width (in eV) 

23 If this parameter is set, then a smoothed variant of the sqrt(E) 

24 expression is used. This parameter sets the energy width of 

25 the smoothing. 

26 The eps parameter is ignored when the width parameter is set. 

27 """ 

28 def __init__(self, 

29 eps: float = 0.05, 

30 width: float = None): 

31 self.width: float | None 

32 if width is not None: 

33 self.width = width / Ha 

34 self.eps = None # Make sure that eps is not used with width 

35 else: 

36 self.width = None 

37 self.eps = eps / Ha 

38 

39 def initialize(self, wfs): 

40 self.wfs = wfs 

41 

42 def initialize_1d(self, ae): 

43 self.ae = ae 

44 

45 def get_description(self): 

46 desc = [] 

47 if self.eps is not None: 

48 desc += ['eps={:.4f} eV'.format(self.eps * Ha)] 

49 if self.width is not None: 

50 desc += ['width={:.4f} eV'.format(self.width * Ha)] 

51 return ', '.join(desc) 

52 

53 def f(self, energy_n): 

54 """Calculate the sqrt(E)-like coefficient. 

55 

56 See the class description for details. 

57 """ 

58 w_n = np.zeros_like(energy_n, dtype=float) 

59 if self.width is None: 

60 flt_n = energy_n > self.eps 

61 w_n[flt_n] = np.sqrt(energy_n[flt_n]) 

62 else: 

63 prefactor = 0.5 * np.sqrt(np.pi * self.width) 

64 rel_energy_n = energy_n / self.width 

65 # Evaluate positive energies 

66 flt_n = energy_n > 0.0 

67 w_n[flt_n] = (np.sqrt(energy_n[flt_n]) 

68 + prefactor * erfcx(np.sqrt(rel_energy_n[flt_n]))) 

69 # Evaluate negative energies 

70 flt_n = np.logical_not(flt_n) 

71 w_n[flt_n] = prefactor * np.exp(rel_energy_n[flt_n]) 

72 return w_n 

73 

74 def get_reference_homo_1d(self): 

75 e_j = np.asarray(self.ae.e_j) 

76 f_j = np.asarray(self.ae.f_j) 

77 homo = np.max(e_j[f_j > 1e-3]) 

78 return homo 

79 

80 def get_reference_lumo_1d(self): 

81 e_j = np.asarray(self.ae.e_j) 

82 f_j = np.asarray(self.ae.f_j) 

83 lumo = np.min(e_j[f_j < 1e-3]) 

84 return lumo 

85 

86 def get_coefficients_1d(self, smooth=False): 

87 homo = self.get_reference_homo_1d() 

88 if smooth: 

89 w_ln = [] 

90 for e_n, f_n in zip(self.ae.e_ln, self.ae.f_ln): 

91 e_n = np.asarray(e_n) 

92 f_n = np.asarray(f_n) 

93 w_n = f_n * K_G * self.f(homo - e_n) 

94 w_ln.append(w_n) 

95 return w_ln 

96 else: 

97 e_j = np.asarray(self.ae.e_j) 

98 f_j = np.asarray(self.ae.f_j) 

99 return f_j * K_G * self.f(homo - e_j) 

100 

101 def get_coefficients_1d_for_lumo_perturbation(self, smooth=False): 

102 homo = self.get_reference_homo_1d() 

103 lumo = self.get_reference_lumo_1d() 

104 e_j = np.asarray(self.ae.e_j) 

105 f_j = np.asarray(self.ae.f_j) 

106 return f_j * K_G * (self.f(lumo - e_j) - self.f(homo - e_j)) 

107 

108 def get_coefficients(self, kpt_u, eref_s): 

109 w_kn = [] 

110 for kpt in kpt_u: 

111 w_n = self.f(eref_s[kpt.s] - kpt.eps_n) 

112 w_n *= kpt.f_n * K_G 

113 w_kn.append(w_n) 

114 return w_kn 

115 

116 def get_coefficients_for_lumo_perturbation(self, kpt_u, eref_s, 

117 eref_lumo_s): 

118 w_kn = [] 

119 for kpt in kpt_u: 

120 w_n = (self.f(eref_lumo_s[kpt.s] - kpt.eps_n) 

121 - self.f(eref_s[kpt.s] - kpt.eps_n)) 

122 w_n *= kpt.f_n * K_G 

123 w_kn.append(w_n) 

124 return w_kn