Coverage for eminus/gth.py: 98.36%

1# SPDX-FileCopyrightText: 2021 The eminus developers

2# SPDX-License-Identifier: Apache-2.0

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

5import math

7import numpy as np

9from . import backend as xp

10from .io import read_gth

11from .utils import Ylm_real

14class GTH:

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

17 Keyword Args:

18 scf: SCF object.

19 """

21 def __init__(self, scf=None):

22 """Initialize the GTH object."""

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

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

25 self.prj2beta = xp.asarray([0]) #: Index matrix to map to the correct projector function.

26 self.betaNL = xp.asarray([0]) #: Atomic-centered projector functions.

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

29 if scf is not None:

30 atoms = scf.atoms

32 # Set up a dictionary with all GTH parameters

33 for ia in range(atoms.Natoms):

34 # Skip if the atom is already in the dictionary

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

36 continue

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

39 # Initialize the non-local potential

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

41 self.NbetaNL = NbetaNL

42 self.prj2beta = prj2beta

43 self.betaNL = betaNL

45 def __getitem__(self, key):

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

47 return self.GTH[key]

49 def __repr__(self):

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

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

54def init_gth_loc(scf, **kwargs):

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

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

59 Args:

60 scf: SCF object.

62 Keyword Args:

63 **kwargs: Throwaway arguments.

65 Returns:

66 Local GTH potential contribution.

67 """

68 atoms = scf.atoms

69 species = set(atoms.atom)

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

72 Vloc = xp.zeros_like(atoms.G2, dtype=complex)

73 for isp in species:

74 psp = scf.gth[isp]

75 rloc = psp["rloc"]

76 Zion = psp["Zion"]

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

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

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

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

82 rlocG2 = atoms.G2 * rloc**2

83 rlocG22 = rlocG2**2

84 exprlocG2 = xp.exp(-0.5 * rlocG2)

85 # Ignore the division by zero for the first elements

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

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

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

89 (2 * math.pi) ** 3

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

91 c1

92 + c2 * (3 - rlocG2)

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

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

95 )

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

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

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

99 )

100

101 # Sum up the structure factor for every species

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

103 for ia in range(atoms.Natoms):

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

105 Sf += atoms.Sf[ia]

106 Vloc += xp.real(atoms.J(Vsp * Sf))

107 return Vloc

108

109

110def init_gth_nonloc(atoms, gth):

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

112

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

114

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

116

117 Args:

118 atoms: Atoms object.

119 gth: GTH object.

120

121 Returns:

122 NbetaNL, prj2beta, and betaNL.

123 """

124 prj2beta = xp.empty((3, atoms.Natoms, 4, 7), dtype=int)

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

126

127 NbetaNL = 0

128 for ia in range(atoms.Natoms):

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

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

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

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

133 NbetaNL += 1

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

135

136 betaNL = []

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

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

139 ibeta = 0

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

141 Gkm = xp.sqrt(atoms.Gk2c[ik])

142 for ia in range(atoms.Natoms):

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

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

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

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

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

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

149 betaNL_ik[:, ibeta] = (

150 (-1j) ** l

151 * Ylm_real(l, m, gk)

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

153 * Sf

154 )

155 ibeta += 1

156 betaNL.append(betaNL_ik)

157 return NbetaNL, prj2beta, betaNL

158

159

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

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

162

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

164

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

166

167 Args:

168 scf: SCF object.

169 ik: k-point index.

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

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

172

173 Returns:

174 Non-local GTH potential contribution.

175 """

176 atoms = scf.atoms

177

178 Vpsi = xp.zeros_like(W[ik][spin], dtype=complex)

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

180 return Vpsi

181

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

183 for ia in range(atoms.Natoms):

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

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

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

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

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

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

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

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

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

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

194 return atoms.O(Vpsi)

195

196

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

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

199

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

201

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

203

204 Args:

205 psp: GTH parameters.

206 l: Angular momentum number.

207 iprj: Nproj_l index.

208 Gm: Magnitude of G-vectors.

209 Omega: Unit cell volume.

210

211 Returns:

212 GTH projector.

213 """

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

215 Gr2 = (Gm * rrl) ** 2

216

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

218 Vprj = prefactor * xp.exp(-0.5 * Gr2)

219

220 if l == 0: # s-channel

221 if iprj == 1:

222 return Vprj

223 if iprj == 2:

224 return 2 / math.sqrt(15) * (3 - Gr2) * Vprj

225 if iprj == 3:

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

227 elif l == 1: # p-channel

228 if iprj == 1:

229 return 1 / math.sqrt(3) * Gm * Vprj

230 if iprj == 2:

231 return 2 / math.sqrt(105) * Gm * (5 - Gr2) * Vprj

232 if iprj == 3:

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

234 elif l == 2: # d-channel

235 if iprj == 1:

236 return 1 / math.sqrt(15) * Gm**2 * Vprj

237 if iprj == 2:

238 return (2 / 3) / math.sqrt(105) * Gm**2 * (7 - Gr2) * Vprj

239 elif l == 3: # f-channel

240 # Only one projector

241 return 1 / math.sqrt(105) * Gm**3 * Vprj

242

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

244 raise ValueError(msg)