Coverage for eminus/potentials.py: 94.20%

1# SPDX-FileCopyrightText: 2021 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

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

4 

5import inspect 

6import math 

7 

8import numpy as np 

9 

10from . import backend as xp 

11from .gth import init_gth_loc 

12from .logger import log 

13 

14 

15def get_pot_defaults(pot): 

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

17 

18 Args: 

19 pot: Potential name. 

20 

21 Returns: 

22 Default parameters and values. 

23 """ 

24 if pot in IMPLEMENTED: 

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

26 return { 

27 param.name: param.default 

28 for param in sig.parameters.values() 

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

30 } 

31 return {} 

32 

33 

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

35 """Harmonic potential. 

36 

37 Can be used for quantum dot calculations. 

38 

39 Args: 

40 scf: SCF object. 

41 

42 Keyword Args: 

43 freq: Harmonic oscillator frequency. 

44 **kwargs: Throwaway arguments. 

45 

46 Returns: 

47 Harmonic potential in real-space. 

48 """ 

49 atoms = scf.atoms 

50 dr = xp.linalg.norm(atoms.r - xp.sum(atoms.a, axis=1) / 2, axis=1) 

51 Vharm = 0.5 * freq**2 * dr**2 

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

53 

54 

55def coulomb(scf, **kwargs): 

56 """All-electron Coulomb potential. 

57 

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

59 

60 Args: 

61 scf: SCF object. 

62 

63 Keyword Args: 

64 **kwargs: Throwaway arguments. 

65 

66 Returns: 

67 Coulomb potential in real-space. 

68 """ 

69 atoms = scf.atoms 

70 

71 Vcoul = xp.zeros_like(atoms.G2, dtype=complex) 

72 for isp in set(atoms.atom): 

73 # Sum up the structure factor for every species 

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

75 Sf = xp.zeros(len(atoms.Sf[0]), dtype=complex) 

76 for ia in range(atoms.Natoms): 

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

78 Sf += atoms.Sf[ia] 

79 Z = atoms.Z[ia] 

80 

81 # Ignore the division by zero for the first elements 

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

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

84 Vsp = -4 * math.pi * Z / atoms.G2 

85 Vsp[0] = 0 

86 Vcoul += xp.real(atoms.J(Vsp * Sf)) 

87 return Vcoul 

88 

89 

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

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

92 

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

94 

95 Args: 

96 scf: SCF object. 

97 

98 Keyword Args: 

99 alpha: Convergence parameter. 

100 **kwargs: Throwaway arguments. 

101 

102 Returns: 

103 Long-range Coulomb potential in real-space. 

104 """ 

105 atoms = scf.atoms 

106 

107 Vcoul = xp.zeros_like(atoms.G2, dtype=complex) 

108 for isp in set(atoms.atom): 

109 # Sum up the structure factor for every species 

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

111 Sf = xp.zeros(len(atoms.Sf[0]), dtype=complex) 

112 for ia in range(atoms.Natoms): 

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

114 Sf += atoms.Sf[ia] 

115 Z = atoms.Z[ia] 

116 

117 # Ignore the division by zero for the first elements 

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

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

120 Vsp = -4 * math.pi * Z * xp.exp(-atoms.G2 / (4 * alpha**2)) / atoms.G2 

121 Vsp[0] = 0 

122 Vcoul += xp.real(atoms.J(Vsp * Sf)) 

123 return Vcoul 

124 

125 

126def ge(scf, **kwargs): 

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

128 

129 Fourier-transformed by Tomas Arias. 

130 

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

132 

133 Args: 

134 scf: SCF object. 

135 

136 Keyword Args: 

137 **kwargs: Throwaway arguments. 

138 

139 Returns: 

140 Germanium pseudopotential in real-space. 

141 """ 

142 atoms = scf.atoms 

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

144 lamda = 18.5 

145 rc = 1.052 

146 Gm = xp.sqrt(atoms.G2) 

147 

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

149 Vps = ( 

150 -2 

151 * math.pi 

152 * xp.exp(-math.pi * Gm / lamda) 

153 * xp.cos(rc * Gm) 

154 * (Gm / lamda) 

155 / (1 - xp.exp(-2 * math.pi * Gm / lamda)) 

156 ) 

157 for n in range(5): 

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

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

160 

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

162 n = xp.arange(1, 5) 

163 Vps[0] = ( 

164 4 

165 * math.pi 

166 * Z 

167 * (1 + math.exp(-lamda * rc)) 

168 * ( 

169 rc**2 / 2 

170 + 1 / lamda**2 * (math.pi**2 / 6 + xp.sum((-1) ** n * xp.exp(-lamda * rc * n) / n**2)) 

171 ) 

172 ) 

173 

174 Sf = xp.sum(atoms.Sf, axis=0) 

175 return atoms.J(Vps * Sf) 

176 

177 

178def init_pot(scf, pot_params=None): 

179 """Handle and initialize potentials. 

180 

181 Args: 

182 scf: SCF object. 

183 

184 Keywords Args: 

185 pot_params: Potential parameters. 

186 

187 Returns: 

188 Potential in real-space. 

189 """ 

190 if pot_params is None: 

191 pot_params = {} 

192 try: 

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

194 except KeyError: 

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

196 raise 

197 return pot 

198 

199 

200#: Map potential names with their respective implementation. 

201IMPLEMENTED = { 

202 "harmonic": harmonic, 

203 "coulomb": coulomb, 

204 "lr": coulomb_lr, 

205 "ge": ge, 

206 "gth": init_gth_loc, 

207}