Coverage for eminus/gth.py: 98.29%

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 # Allow creating an empty instance (used when loading GTH objects from JSON files) 

21 if scf is not None: 

22 atoms = scf.atoms 

23 

24 # Set up a dictionary with all GTH parameters 

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

26 for ia in range(atoms.Natoms): 

27 # Skip if the atom is already in the dictionary 

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

29 continue 

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

31 

32 # Initialize the non-local potential 

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

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

35 self.prj2beta = prj2beta #: Index matrix to map to the correct projector function. 

36 self.betaNL = betaNL #: Atomic-centered projector functions. 

37 

38 def __getitem__(self, key): 

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

40 return self.GTH[key] 

41 

42 def __repr__(self): 

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

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

45 

46 

47def init_gth_loc(scf, **kwargs): 

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

49 

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

51 

52 Args: 

53 scf: SCF object. 

54 

55 Keyword Args: 

56 **kwargs: Throwaway arguments. 

57 

58 Returns: 

59 Local GTH potential contribution. 

60 """ 

61 atoms = scf.atoms 

62 species = set(atoms.atom) 

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

64 

65 Vloc = np.zeros_like(atoms.G2) 

66 for isp in species: 

67 psp = scf.gth[isp] 

68 rloc = psp["rloc"] 

69 Zion = psp["Zion"] 

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

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

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

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

74 

75 rlocG2 = atoms.G2 * rloc**2 

76 rlocG22 = rlocG2**2 

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

78 # Ignore the division by zero for the first elements 

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

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

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

82 (2 * np.pi) ** 3 

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

84 c1 

85 + c2 * (3 - rlocG2) 

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

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

88 ) 

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

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

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

92 ) 

93 

94 # Sum up the structure factor for every species 

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

96 for ia in range(atoms.Natoms): 

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

98 Sf += atoms.Sf[ia] 

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

100 return Vloc 

101 

102 

103def init_gth_nonloc(atoms, gth): 

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

105 

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

107 

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

109 

110 Args: 

111 atoms: Atoms object. 

112 gth: GTH object. 

113 

114 Returns: 

115 NbetaNL, prj2beta, and betaNL. 

116 """ 

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

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

119 

120 NbetaNL = 0 

121 for ia in range(atoms.Natoms): 

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

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

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

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

126 NbetaNL += 1 

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

128 

129 betaNL = [] 

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

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

132 ibeta = 0 

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

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

135 for ia in range(atoms.Natoms): 

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

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

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

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

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

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

142 betaNL_ik[:, ibeta] = ( 

143 (-1j) ** l 

144 * Ylm_real(l, m, gk) 

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

146 * Sf 

147 ) 

148 ibeta += 1 

149 betaNL.append(betaNL_ik) 

150 return NbetaNL, prj2beta, betaNL 

151 

152 

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

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

155 

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

157 

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

159 

160 Args: 

161 scf: SCF object. 

162 ik: k-point index. 

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

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

165 

166 Returns: 

167 Non-local GTH potential contribution. 

168 """ 

169 atoms = scf.atoms 

170 

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

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

173 return Vpsi 

174 

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

176 for ia in range(atoms.Natoms): 

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

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

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

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

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

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

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

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

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

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

187 return atoms.O(Vpsi) 

188 

189 

190def eval_proj_G(psp, l, iprj, Gm, Omega): # noqa: PLR0911 

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

192 

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

194 

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

196 

197 Args: 

198 psp: GTH parameters. 

199 l: Angular momentum number. 

200 iprj: Nproj_l index. 

201 Gm: Magnitude of G-vectors. 

202 Omega: Unit cell volume. 

203 

204 Returns: 

205 GTH projector. 

206 """ 

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

208 Gr2 = (Gm * rrl) ** 2 

209 

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

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

212 

213 if l == 0: # s-channel 

214 if iprj == 1: 

215 return Vprj 

216 if iprj == 2: 

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

218 if iprj == 3: 

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

220 elif l == 1: # p-channel 

221 if iprj == 1: 

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

223 if iprj == 2: 

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

225 if iprj == 3: 

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

227 elif l == 2: # d-channel 

228 if iprj == 1: 

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

230 if iprj == 2: 

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

232 elif l == 3: # f-channel 

233 # Only one projector 

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

235 

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

237 raise ValueError(msg)