Coverage for eminus/gth.py: 98.33%

1# SPDX-FileCopyrightText: 2021 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""Utilities to use Goedecker, Teter, and Hutter pseudopotentials.""" 

4 

5import numpy as np 

6 

7from .io import read_gth 

8from .utils import Ylm_real 

9 

10 

11class GTH: 

12 """GTH object that holds non-local projectors and parameters for all atoms in the SCF object. 

13 

14 Keyword Args: 

15 scf: SCF object. 

16 """ 

17 

18 def __init__(self, scf=None): 

19 """Initialize the GTH object.""" 

20 self.GTH = {} #: Dictionary with GTH parameters per atom type. 

21 self.NbetaNL = 0 #: Number of projector functions for the non-local potential. 

22 self.prj2beta = np.array([0]) #: Index matrix to map to the correct projector function. 

23 self.betaNL = np.array([0]) #: Atomic-centered projector functions. 

24 

25 # Allow creating an empty instance (used when loading GTH objects from JSON files) 

26 if scf is not None: 

27 atoms = scf.atoms 

28 

29 # Set up a dictionary with all GTH parameters 

30 for ia in range(atoms.Natoms): 

31 # Skip if the atom is already in the dictionary 

32 if atoms.atom[ia] in self.GTH: 

33 continue 

34 self.GTH[atoms.atom[ia]] = read_gth(atoms.atom[ia], atoms.Z[ia], psp_path=scf.psp) 

35 

36 # Initialize the non-local potential 

37 NbetaNL, prj2beta, betaNL = init_gth_nonloc(atoms, self) 

38 self.NbetaNL = NbetaNL 

39 self.prj2beta = prj2beta 

40 self.betaNL = betaNL 

41 

42 def __getitem__(self, key): 

43 """Allow accessing the GTH parameters of an atom by indexing the GTH object.""" 

44 return self.GTH[key] 

45 

46 def __repr__(self): 

47 """Print a short overview over the values stored in the GTH object.""" 

48 return f"NbetaNL: {self.NbetaNL}\nGTH values for: {', '.join(list(self.GTH))}" 

49 

50 

51def init_gth_loc(scf, **kwargs): 

52 """Initialize parameters to calculate local contributions of GTH pseudopotentials. 

53 

54 Reference: Phys. Rev. B 54, 1703. 

55 

56 Args: 

57 scf: SCF object. 

58 

59 Keyword Args: 

60 **kwargs: Throwaway arguments. 

61 

62 Returns: 

63 Local GTH potential contribution. 

64 """ 

65 atoms = scf.atoms 

66 species = set(atoms.atom) 

67 omega = 1 # Normally this would be det(atoms.a), but Arias notation is off by this factor 

68 

69 Vloc = np.zeros_like(atoms.G2) 

70 for isp in species: 

71 psp = scf.gth[isp] 

72 rloc = psp["rloc"] 

73 Zion = psp["Zion"] 

74 c1 = psp["cloc"][0] 

75 c2 = psp["cloc"][1] 

76 c3 = psp["cloc"][2] 

77 c4 = psp["cloc"][3] 

78 

79 rlocG2 = atoms.G2 * rloc**2 

80 rlocG22 = rlocG2**2 

81 exprlocG2 = np.exp(-0.5 * rlocG2) 

82 # Ignore the division by zero for the first elements 

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

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

85 Vsp = -4 * np.pi * Zion / omega * exprlocG2 / atoms.G2 + np.sqrt( 

86 (2 * np.pi) ** 3 

87 ) * rloc**3 / omega * exprlocG2 * ( 

88 c1 

89 + c2 * (3 - rlocG2) 

90 + c3 * (15 - 10 * rlocG2 + rlocG22) 

91 + c4 * (105 - 105 * rlocG2 + 21 * rlocG22 - rlocG2**3) 

92 ) 

93 # Special case for G=(0,0,0), same as in QE 

94 Vsp[0] = 2 * np.pi * Zion * rloc**2 + (2 * np.pi) ** 1.5 * rloc**3 * ( 

95 c1 + 3 * c2 + 15 * c3 + 105 * c4 

96 ) 

97 

98 # Sum up the structure factor for every species 

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

100 for ia in range(atoms.Natoms): 

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

102 Sf += atoms.Sf[ia] 

103 Vloc += np.real(atoms.J(Vsp * Sf)) 

104 return Vloc 

105 

106 

107def init_gth_nonloc(atoms, gth): 

108 """Initialize parameters to calculate non-local contributions of GTH pseudopotentials. 

109 

110 Adapted from https://github.com/f-fathurrahman/PWDFT.jl/blob/master/src/PsPotNL.jl 

111 

112 Reference: Phys. Rev. B 54, 1703. 

113 

114 Args: 

115 atoms: Atoms object. 

116 gth: GTH object. 

117 

118 Returns: 

119 NbetaNL, prj2beta, and betaNL. 

120 """ 

121 prj2beta = np.empty((3, atoms.Natoms, 4, 7), dtype=int) 

122 prj2beta[:] = -1 # Set to an invalid index 

123 

124 NbetaNL = 0 

125 for ia in range(atoms.Natoms): 

126 psp = gth[atoms.atom[ia]] 

127 for l in range(psp["lmax"]): 

128 for m in range(-l, l + 1): 

129 for iprj in range(psp["Nproj_l"][l]): 

130 NbetaNL += 1 

131 prj2beta[iprj, ia, l, m + psp["lmax"] - 1] = NbetaNL 

132 

133 betaNL = [] 

134 for ik in range(atoms.kpts.Nk): 

135 betaNL_ik = np.empty((len(atoms.Gk2c[ik]), NbetaNL), dtype=complex) 

136 ibeta = 0 

137 gk = atoms.G[atoms.active[ik]] + atoms.kpts.k[ik] 

138 Gkm = np.sqrt(atoms.Gk2c[ik]) 

139 for ia in range(atoms.Natoms): 

140 # It is important to transform the structure factor to make both notations compatible 

141 Sf = atoms.Idag(atoms.J(atoms.Sf[ia], ik), ik) 

142 psp = gth[atoms.atom[ia]] 

143 for l in range(psp["lmax"]): 

144 for m in range(-l, l + 1): 

145 for iprj in range(psp["Nproj_l"][l]): 

146 betaNL_ik[:, ibeta] = ( 

147 (-1j) ** l 

148 * Ylm_real(l, m, gk) 

149 * eval_proj_G(psp, l, iprj + 1, Gkm, atoms.Omega) 

150 * Sf 

151 ) 

152 ibeta += 1 

153 betaNL.append(betaNL_ik) 

154 return NbetaNL, prj2beta, betaNL 

155 

156 

157def calc_Vnonloc(scf, ik, spin, W): 

158 """Calculate the non-local pseudopotential, applied on the basis functions W. 

159 

160 Adapted from https://github.com/f-fathurrahman/PWDFT.jl/blob/master/src/op_V_Ps_nloc.jl 

161 

162 Reference: Phys. Rev. B 54, 1703. 

163 

164 Args: 

165 scf: SCF object. 

166 ik: k-point index. 

167 spin: Spin variable to track whether to do the calculation for spin up or down. 

168 W: Expansion coefficients of unconstrained wave functions in reciprocal space. 

169 

170 Returns: 

171 Non-local GTH potential contribution. 

172 """ 

173 atoms = scf.atoms 

174 

175 Vpsi = np.zeros_like(W[ik][spin], dtype=complex) 

176 if scf.pot != "gth" or scf.gth.NbetaNL == 0: # Only calculate the non-local part if necessary 

177 return Vpsi 

178 

179 betaNL_psi = (W[ik][spin].conj().T @ scf.gth.betaNL[ik]).conj() 

180 for ia in range(atoms.Natoms): 

181 psp = scf.gth[atoms.atom[ia]] 

182 for l in range(psp["lmax"]): 

183 for m in range(-l, l + 1): 

184 for iprj in range(psp["Nproj_l"][l]): 

185 ibeta = scf.gth.prj2beta[iprj, ia, l, m + psp["lmax"] - 1] - 1 

186 for jprj in range(psp["Nproj_l"][l]): 

187 jbeta = scf.gth.prj2beta[jprj, ia, l, m + psp["lmax"] - 1] - 1 

188 hij = psp["h"][l, iprj, jprj] 

189 Vpsi += hij * betaNL_psi[:, jbeta] * scf.gth.betaNL[ik][:, ibeta, None] 

190 # We have to multiply with the cell volume, because of different orthogonalization methods 

191 return atoms.O(Vpsi) 

192 

193 

194def eval_proj_G(psp, l, iprj, Gm, Omega): 

195 """Evaluate GTH projector functions in G-space. 

196 

197 Adapted from https://github.com/f-fathurrahman/PWDFT.jl/blob/master/src/PsPot_GTH.jl 

198 

199 Reference: Phys. Rev. B 54, 1703. 

200 

201 Args: 

202 psp: GTH parameters. 

203 l: Angular momentum number. 

204 iprj: Nproj_l index. 

205 Gm: Magnitude of G-vectors. 

206 Omega: Unit cell volume. 

207 

208 Returns: 

209 GTH projector. 

210 """ 

211 rrl = psp["rp"][l] 

212 Gr2 = (Gm * rrl) ** 2 

213 

214 prefactor = 4 * np.pi ** (5 / 4) * np.sqrt(2 ** (l + 1) * rrl ** (2 * l + 3) / Omega) 

215 Vprj = prefactor * np.exp(-0.5 * Gr2) 

216 

217 if l == 0: # s-channel 

218 if iprj == 1: 

219 return Vprj 

220 if iprj == 2: 

221 return 2 / np.sqrt(15) * (3 - Gr2) * Vprj 

222 if iprj == 3: 

223 return (4 / 3) / np.sqrt(105) * (15 - 10 * Gr2 + Gr2**2) * Vprj 

224 elif l == 1: # p-channel 

225 if iprj == 1: 

226 return 1 / np.sqrt(3) * Gm * Vprj 

227 if iprj == 2: 

228 return 2 / np.sqrt(105) * Gm * (5 - Gr2) * Vprj 

229 if iprj == 3: 

230 return (4 / 3) / np.sqrt(1155) * Gm * (35 - 14 * Gr2 + Gr2**2) * Vprj 

231 elif l == 2: # d-channel 

232 if iprj == 1: 

233 return 1 / np.sqrt(15) * Gm**2 * Vprj 

234 if iprj == 2: 

235 return (2 / 3) / np.sqrt(105) * Gm**2 * (7 - Gr2) * Vprj 

236 elif l == 3: # f-channel 

237 # Only one projector 

238 return 1 / np.sqrt(105) * Gm**3 * Vprj 

239 

240 msg = f"No projector found for l={l}." 

241 raise ValueError(msg)