Coverage for eminus/potentials.py: 94.12%

1# SPDX-FileCopyrightText: 2021 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""Collection of miscellaneous potentials.""" 

4 

5import inspect 

6 

7import numpy as np 

8from scipy.linalg import norm 

9 

10from .gth import init_gth_loc 

11from .logger import log 

12 

13 

14def get_pot_defaults(pot): 

15 """Get the default parameters and values for a given potential. 

16 

17 Args: 

18 pot: Potential name. 

19 

20 Returns: 

21 Default parameters and values. 

22 """ 

23 if pot in IMPLEMENTED: 

24 sig = inspect.signature(IMPLEMENTED[pot]) 

25 return { 

26 param.name: param.default 

27 for param in sig.parameters.values() 

28 if param.default is not inspect.Parameter.empty 

29 } 

30 return {} 

31 

32 

33def harmonic(scf, freq=2, **kwargs): 

34 """Harmonic potential. 

35 

36 Can be used for quantum dot calculations. 

37 

38 Args: 

39 scf: SCF object. 

40 

41 Keyword Args: 

42 freq: Harmonic oscillator frequency. 

43 **kwargs: Throwaway arguments. 

44 

45 Returns: 

46 Harmonic potential in real-space. 

47 """ 

48 atoms = scf.atoms 

49 dr = norm(atoms.r - np.sum(atoms.a, axis=1) / 2, axis=1) 

50 Vharm = 0.5 * freq**2 * dr**2 

51 return atoms.Jdag(atoms.O(atoms.J(Vharm))) 

52 

53 

54def coulomb(scf, **kwargs): 

55 """All-electron Coulomb potential. 

56 

57 Reference: Bull. Lebedev Phys. Inst. 42, 329. 

58 

59 Args: 

60 scf: SCF object. 

61 

62 Keyword Args: 

63 **kwargs: Throwaway arguments. 

64 

65 Returns: 

66 Coulomb potential in real-space. 

67 """ 

68 atoms = scf.atoms 

69 

70 Vcoul = np.zeros_like(atoms.G2) 

71 for isp in set(atoms.atom): 

72 # Sum up the structure factor for every species 

73 # Also get the charge, assuming all species have the same charge 

74 Sf = np.zeros(len(atoms.Sf[0]), dtype=complex) 

75 for ia in range(atoms.Natoms): 

76 if atoms.atom[ia] == isp: 

77 Sf += atoms.Sf[ia] 

78 Z = atoms.Z[ia] 

79 

80 # Ignore the division by zero for the first elements 

81 # One could do some proper indexing with [1:] but indexing is slow 

82 with np.errstate(divide="ignore", invalid="ignore"): 

83 Vsp = -4 * np.pi * Z / atoms.G2 

84 Vsp[0] = 0 

85 Vcoul += np.real(atoms.J(Vsp * Sf)) 

86 return Vcoul 

87 

88 

89def coulomb_lr(scf, alpha=100, **kwargs): 

90 """Long-range all-electron Coulomb potential. 

91 

92 Reference: J. Comput. Phys. 117, 171. 

93 

94 Args: 

95 scf: SCF object. 

96 

97 Keyword Args: 

98 alpha: Convergence parameter. 

99 **kwargs: Throwaway arguments. 

100 

101 Returns: 

102 Long-range Coulomb potential in real-space. 

103 """ 

104 atoms = scf.atoms 

105 

106 Vcoul = np.zeros_like(atoms.G2) 

107 for isp in set(atoms.atom): 

108 # Sum up the structure factor for every species 

109 # Also get the charge, assuming all species have the same charge 

110 Sf = np.zeros(len(atoms.Sf[0]), dtype=complex) 

111 for ia in range(atoms.Natoms): 

112 if atoms.atom[ia] == isp: 

113 Sf += atoms.Sf[ia] 

114 Z = atoms.Z[ia] 

115 

116 # Ignore the division by zero for the first elements 

117 # One could do some proper indexing with [1:] but indexing is slow 

118 with np.errstate(divide="ignore", invalid="ignore"): 

119 Vsp = -4 * np.pi * Z * np.exp(-atoms.G2 / (4 * alpha**2)) / atoms.G2 

120 Vsp[0] = 0 

121 Vcoul += np.real(atoms.J(Vsp * Sf)) 

122 return Vcoul 

123 

124 

125def ge(scf, **kwargs): 

126 """Starkloff-Joannopoulos local pseudopotential for germanium. 

127 

128 Fourier-transformed by Tomas Arias. 

129 

130 Reference: Phys. Rev. B 16, 5212. 

131 

132 Args: 

133 scf: SCF object. 

134 

135 Keyword Args: 

136 **kwargs: Throwaway arguments. 

137 

138 Returns: 

139 Germanium pseudopotential in real-space. 

140 """ 

141 atoms = scf.atoms 

142 Z = 4 # This potential should only be used for germanium 

143 lamda = 18.5 

144 rc = 1.052 

145 Gm = np.sqrt(atoms.G2) 

146 

147 with np.errstate(divide="ignore", invalid="ignore"): 

148 Vps = ( 

149 -2 

150 * np.pi 

151 * np.exp(-np.pi * Gm / lamda) 

152 * np.cos(rc * Gm) 

153 * (Gm / lamda) 

154 / (1 - np.exp(-2 * np.pi * Gm / lamda)) 

155 ) 

156 for n in range(5): 

157 Vps = Vps + (-1) ** n * np.exp(-lamda * rc * n) / (1 + (n * lamda / Gm) ** 2) 

158 Vps = Vps * 4 * np.pi * Z / Gm**2 * (1 + np.exp(-lamda * rc)) - 4 * np.pi * Z / Gm**2 

159 

160 # Special case for G=(0,0,0) 

161 n = np.arange(1, 5) 

162 Vps[0] = ( 

163 4 

164 * np.pi 

165 * Z 

166 * (1 + np.exp(-lamda * rc)) 

167 * ( 

168 rc**2 / 2 

169 + 1 / lamda**2 * (np.pi**2 / 6 + np.sum((-1) ** n * np.exp(-lamda * rc * n) / n**2)) 

170 ) 

171 ) 

172 

173 Sf = np.sum(atoms.Sf, axis=0) 

174 return atoms.J(Vps * Sf) 

175 

176 

177def init_pot(scf, pot_params=None): 

178 """Handle and initialize potentials. 

179 

180 Args: 

181 scf: SCF object. 

182 

183 Keywords Args: 

184 pot_params: Potential parameters. 

185 

186 Returns: 

187 Potential in real-space. 

188 """ 

189 if pot_params is None: 

190 pot_params = {} 

191 try: 

192 pot = IMPLEMENTED[scf.pot](scf, **pot_params) 

193 except KeyError: 

194 log.exception(f'No potential found for "{scf.pot}".') 

195 raise 

196 return pot 

197 

198 

199#: Map potential names with their respective implementation. 

200IMPLEMENTED = { 

201 "harmonic": harmonic, 

202 "coulomb": coulomb, 

203 "lr": coulomb_lr, 

204 "ge": ge, 

205 "gth": init_gth_loc, 

206}