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): 

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 Returns: 

56 Local GTH potential contribution. 

57 """ 

58 atoms = scf.atoms 

59 species = set(atoms.atom) 

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

61 

62 Vloc = np.zeros_like(atoms.G2) 

63 for isp in species: 

64 psp = scf.gth[isp] 

65 rloc = psp["rloc"] 

66 Zion = psp["Zion"] 

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

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

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

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

71 

72 rlocG2 = atoms.G2 * rloc**2 

73 rlocG22 = rlocG2**2 

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

75 # Ignore the division by zero for the first elements 

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

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

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

79 (2 * np.pi) ** 3 

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

81 c1 

82 + c2 * (3 - rlocG2) 

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

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

85 ) 

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

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

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

89 ) 

90 

91 # Sum up the structure factor for every species 

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

93 for ia in range(atoms.Natoms): 

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

95 Sf += atoms.Sf[ia] 

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

97 return Vloc 

98 

99 

100def init_gth_nonloc(atoms, gth): 

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

102 

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

104 

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

106 

107 Args: 

108 atoms: Atoms object. 

109 gth: GTH object. 

110 

111 Returns: 

112 NbetaNL, prj2beta, and betaNL. 

113 """ 

114 prj2beta = np.zeros((3, atoms.Natoms, 4, 7), dtype=int) 

115 prj2beta += -1 # Set to an invalid index 

116 

117 NbetaNL = 0 

118 for ia in range(atoms.Natoms): 

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

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

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

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

123 NbetaNL += 1 

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

125 

126 betaNL = [] 

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

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

129 ibeta = 0 

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

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

132 for ia in range(atoms.Natoms): 

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

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

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

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

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

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

139 betaNL_ik[:, ibeta] = ( 

140 (-1j) ** l 

141 * Ylm_real(l, m, gk) 

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

143 * Sf 

144 ) 

145 ibeta += 1 

146 betaNL.append(betaNL_ik) 

147 return NbetaNL, prj2beta, betaNL 

148 

149 

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

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

152 

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

154 

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

156 

157 Args: 

158 scf: SCF object. 

159 ik: k-point index. 

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

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

162 

163 Returns: 

164 Non-local GTH potential contribution. 

165 """ 

166 atoms = scf.atoms 

167 

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

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

170 return Vpsi 

171 

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

173 for ia in range(atoms.Natoms): 

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

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

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

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

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

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

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

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

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

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

184 return atoms.O(Vpsi) 

185 

186 

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

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

189 

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

191 

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

193 

194 Args: 

195 psp: GTH parameters. 

196 l: Angular momentum number. 

197 iprj: Nproj_l index. 

198 Gm: Magnitude of G-vectors. 

199 Omega: Unit cell volume. 

200 

201 Returns: 

202 GTH projector. 

203 """ 

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

205 Gr2 = (Gm * rrl) ** 2 

206 

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

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

209 

210 if l == 0: # s-channel 

211 if iprj == 1: 

212 return Vprj 

213 if iprj == 2: 

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

215 if iprj == 3: 

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

217 elif l == 1: # p-channel 

218 if iprj == 1: 

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

220 if iprj == 2: 

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

222 if iprj == 3: 

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

224 elif l == 2: # d-channel 

225 if iprj == 1: 

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

227 if iprj == 2: 

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

229 elif l == 3: # f-channel 

230 # Only one projector 

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

232 

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

234 raise ValueError(msg)