Coverage for gpaw/hgh.py: 79%

1import hashlib

2from typing import Dict

4import numpy as np

5from ase.data import atomic_numbers

7from gpaw.utilities import pack_hermitian

8from gpaw.atom.radialgd import AERadialGridDescriptor

9from gpaw.atom.configurations import configurations

10from gpaw.pseudopotential import PseudoPotential, get_radial_hartree_energy

13setups: Dict[str, 'HGHParameterSet'] = {} # Filled out during parsing below

14sc_setups: Dict[str, 'HGHParameterSet'] = {} # Semicore

17# Tabulated values of Gamma(m + 1/2)

18half_integer_gamma = [np.sqrt(np.pi)]

19for m in range(20):

20 half_integer_gamma.append(half_integer_gamma[m] * (m + 0.5))

23class HGHSetupData:

24 """Setup-compatible class implementing HGH pseudopotential.

26 To the PAW code this will appear as a legit PAW setup, but is

27 in fact considerably simpler. In particular, all-electron and

28 pseudo partial waves are all zero, and compensation charges do not

29 depend on environment.

31 A HGH setup has the following form::

33 ----

34 \

35 V = Vlocal + ) | p > h < p |

36 / i ij j

37 ----

38 ij

40 Vlocal contains a short-range term which is Gaussian-shaped and

41 implemented as vbar of a PAW setup, along with a long-range term

42 which goes like 1/r and is implemented in terms of a compensation

43 charge.

45 The non-local part contains KB projector functions which are

46 essentially similar to those in PAW, while h_ij are constants.

47 h_ij are provided by setting the K_p variable of the normal

48 setup.

50 Most other properties of a PAW setup do not exist for HGH setups, for

51 which reason they are generally set to zero:

53 * All-electron partial waves: always zero

54 * Pseudo partial waves: always zero

55 * Projectors: HGH projectors

56 * Zero potential (vbar): Gaussian times polynomial

57 * Compensation charges: One Gaussian-shaped spherically symmetric charge

58 * All-electron core density: Delta function corresponding to core electron

59 charge

60 * Pseudo core density: always zero

61 * Pseudo valence density: always zero

62 * PS/AE Kinetic energy density: always zero

63 * The mysterious constants K_p of a setup correspond to h_ij.

65 Note that since the pseudo partial waves are set to zero,

66 initialization of atomic orbitals requires loading a custom basis

67 set.

69 Absolute energies become numerically large since no atomic

70 reference is subtracted.

71 """

72 def __init__(self, hghdata):

73 if isinstance(hghdata, str):

74 symbol = hghdata

75 if symbol.endswith('.sc'):

76 hghdata = sc_setups[symbol[:-3]]

77 else:

78 hghdata = setups[symbol]

79 self.hghdata = hghdata

81 chemsymbol = hghdata.symbol

82 if '.' in chemsymbol:

83 chemsymbol, sc = chemsymbol.split('.')

84 assert sc == 'sc'

85 self.symbol = chemsymbol

86 self.type = hghdata.symbol

87 self.name = 'LDA'

88 self.initialize_setup_data()

90 def initialize_setup_data(self):

91 hghdata = self.hghdata

92 beta = 0.1

93 N = 450

94 rgd = AERadialGridDescriptor(beta / N, 1.0 / N, N,

95 default_spline_points=100)

96 # rgd = EquidistantRadialGridDescriptor(0.001, 10000)

97 self.rgd = rgd

99 self.Z = hghdata.Z

100 self.Nc = hghdata.Z - hghdata.Nv

101 self.Nv = hghdata.Nv

102 self.rcgauss = np.sqrt(2.0) * hghdata.rloc

103

104 threshold = 1e-8

105 if len(hghdata.c_n) > 0:

106 vloc_g = create_local_shortrange_potential(rgd.r_g, hghdata.rloc,

107 hghdata.c_n)

108 gcutvbar, rcutvbar = self.find_cutoff(rgd.r_g, rgd.dr_g, vloc_g,

109 threshold)

110 self.vbar_g = np.sqrt(4.0 * np.pi) * vloc_g[:gcutvbar]

111 else:

112 rcutvbar = 0.5

113 gcutvbar = rgd.ceil(rcutvbar)

114 self.vbar_g = np.zeros(gcutvbar)

115

116 nj = sum([v.nn for v in hghdata.v_l])

117 if nj == 0:

118 nj = 1 # Code assumes nj > 0 elsewhere, we fill out with zeroes

119

120 if not hghdata.v_l:

121 # No projectors. But the remaining code assumes that everything

122 # has projectors! We'll just add the zero function then

123 hghdata.v_l = [VNonLocal(0, 0.01, [0.])]

124

125 n_j = []

126 l_j = []

127

128 # j ordering is significant, must be nl rather than ln

129 for n, l in self.hghdata.nl_iter():

130 n_j.append(n + 1) # Note: actual n must be positive!

131 l_j.append(l)

132 assert nj == len(n_j)

133 self.nj = nj

134 self.l_j = l_j

135 self.l_orb_J = l_j

136 self.n_j = n_j

137

138 self.rcut_j = []

139 self.pt_jg = []

140

141 for n, l in zip(n_j, l_j):

142 # Note: even pseudopotentials without projectors will get one

143 # projector, but the coefficients h_ij should be zero so it

144 # doesn't matter

145 pt_g = create_hgh_projector(rgd.r_g, l, n, hghdata.v_l[l].r0)

146 norm = np.sqrt(np.dot(rgd.dr_g, pt_g**2 * rgd.r_g**2))

147 assert np.abs(1 - norm) < 1e-5, str(1 - norm)

148 gcut, rcut = self.find_cutoff(rgd.r_g, rgd.dr_g, pt_g, threshold)

149 if rcut < 0.5:

150 rcut = 0.5

151 gcut = rgd.ceil(rcut)

152 pt_g = pt_g[:gcut].copy()

153 rcut = max(rcut, 0.5)

154 self.rcut_j.append(rcut)

155 self.pt_jg.append(pt_g)

156

157 # This is the correct magnitude of the otherwise normalized

158 # compensation charge

159 self.Delta0 = -self.Nv / np.sqrt(4.0 * np.pi)

160

161 f_ln = self.hghdata.get_occupation_numbers()

162 f_j = [0] * nj

163 for j, (n, l) in enumerate(self.hghdata.nl_iter()):

164 try:

165 f_j[j] = f_ln[l][n]

166 except IndexError:

167 pass

168 self.f_ln = f_ln

169 self.f_j = f_j

170

171 r_g, lcomp, ghat = self.get_compensation_charge_functions()

172 assert lcomp == [0] and len(ghat) == 1

173 renormalized_ghat = self.Nv / (4.0 * np.pi) * ghat[0]

174 self.Eh_compcharge = get_radial_hartree_energy(r_g, renormalized_ghat)

175

176 def find_cutoff(self, r_g, dr_g, f_g, sqrtailnorm=1e-5):

177 g = len(r_g)

178 acc_sqrnorm = 0.0

179 while acc_sqrnorm <= sqrtailnorm:

180 g -= 1

181 acc_sqrnorm += (r_g[g] * f_g[g])**2.0 * dr_g[g]

182 if r_g[g] < 0.5: # XXX

183 return g, r_g[g]

184 return g, r_g[g]

185

186 def expand_hamiltonian_matrix(self):

187 """Construct K_p from individual h_nn for each l."""

188 ni = sum([2 * l + 1 for l in self.l_j])

189

190 H_ii = np.zeros((ni, ni))

191

192 # The H_ii used in gpaw is much larger and more general than the one

193 # required for HGH pseudopotentials. This means a lot of the elements

194 # must be assigned the same value. Not a performance issue though,

195 # since these are small matrices

196 M1start = 0

197 for n1, l1 in self.hghdata.nl_iter():

198 M1end = M1start + 2 * l1 + 1

199 M2start = 0

200 v = self.hghdata.v_l[l1]

201 for n2, l2 in self.hghdata.nl_iter():

202 M2end = M2start + 2 * l2 + 1

203 if l1 == l2:

204 h_nn = v.expand_hamiltonian_diagonal()

205 H_mm = np.identity(M2end - M2start) * h_nn[n1, n2]

206 H_ii[M1start:M1end, M2start:M2end] += H_mm

207 M2start = M2end

208 M1start = M1end

209 K_p = pack_hermitian(H_ii)

210 return K_p

211

212 def __str__(self):

213 return "HGHSetupData('%s')" % self.type

214

215 def __repr__(self):

216 return self.__str__()

217

218 def print_info(self, text, _setup):

219 self.hghdata.print_info(text)

220

221 def plot(self):

222 """Plot localized functions of HGH setup."""

223 import matplotlib.pyplot as plt

224 rgd = self.rgd

225

226 plt.subplot(211) # vbar, compensation charge

227 gcutvbar = len(self.vbar_g)

228 plt.plot(rgd.r_g[:gcutvbar], self.vbar_g, 'r', label='vloc',

229 linewidth=3)

230 rcc, gcc = self.get_compensation_charge_functions()

231 gcc = gcc[0]

232

233 plt.plot(rcc, gcc * self.Delta0, 'b--',

234 label='Comp charge [arb. unit]',

235 linewidth=3)

236 plt.legend(loc='best')

237

238 plt.subplot(212) # projectors

239 for j, (n, l, pt_g) in enumerate(zip(self.n_j, self.l_j, self.pt_jg)):

240 label = 'n=%d, l=%d' % (n, l)

241 plt.ylabel('$p_n^l(r)$')

242 ng = len(pt_g)

243 r_g = rgd.r_g[:ng]

244 plt.plot(r_g, pt_g, label=label)

245 plt.legend()

246

247 def create_basis_functions(self):

248 from gpaw.pseudopotential import generate_basis_functions

249 return generate_basis_functions(self)

250

251 def get_compensation_charge_functions(self):

252 alpha = self.rcgauss**-2

253

254 rcutgauss = self.rcgauss * 5.0

255 # smaller values break charge conservation

256

257 r = np.linspace(0.0, rcutgauss, 100)

258 g = alpha**1.5 * np.exp(-alpha * r**2) * 4.0 / np.sqrt(np.pi)

259 g[-1] = 0.0

260 return r, [0], [g]

261

262 def get_projectors(self):

263 # XXX equal-range projectors still required for some reason

264 maxlen = max([len(pt_g) for pt_g in self.pt_jg])

265 pt_j = []

266 for l, pt1_g in zip(self.l_j, self.pt_jg):

267 pt2_g = self.rgd.zeros()[:maxlen]

268 pt2_g[:len(pt1_g)] = pt1_g

269 pt_j.append(self.rgd.spline(pt2_g, self.rgd.r_g[maxlen - 1], l))

270 return pt_j

271

272 def get_local_potential(self):

273 n = len(self.vbar_g)

274 return self.rgd.spline(self.vbar_g, self.rgd.r_g[n - 1], l=0)

275

276 def build(self, xcfunc, lmax, basis, filter=None):

277 if basis is None:

278 basis = self.create_basis_functions()

279 setup = PseudoPotential(self, basis)

280 setup.fingerprint = hashlib.md5(str(self.hghdata).encode()).hexdigest()

281 return setup

282

283

284def create_local_shortrange_potential(r_g, rloc, c_n):

285 rr_g = r_g / rloc # "Relative r"

286 rr2_g = rr_g**2

287 rr4_g = rr2_g**2

288 rr6_g = rr4_g * rr2_g

289

290 gaussianpart = np.exp(-.5 * rr2_g)

291 polypart = np.zeros(r_g.shape)

292 for c, rrn_g in zip(c_n, [1, rr2_g, rr4_g, rr6_g]):

293 polypart += c * rrn_g

294

295 vloc_g = gaussianpart * polypart

296 return vloc_g

297

298

299def create_hgh_projector(r_g, l, n, r0):

300 poly_g = r_g**(l + 2 * (n - 1))

301 gauss_g = np.exp(-.5 * r_g**2 / r0**2)

302 A = r0**(l + (4 * n - 1) / 2.0)

303 assert (4 * n - 1) % 2 == 1

304 B = half_integer_gamma[l + (4 * n - 1) // 2]**.5

305 pt_g = 2.**.5 / A / B * poly_g * gauss_g

306 return pt_g

307

308

309# Coefficients determining off-diagonal elements of h_nn for l = 0...2

310# given the diagonal elements

311hcoefs_l = [

312 [-.5 * (3. / 5.)**.5, .5 * (5. / 21.)**.5, -.5 * (100. / 63.)**.5],

313 [-.5 * (5. / 7.)**.5, 1. / 6. * (35. / 11.)**.5, -1. / 6. * 14. / 11.**.5],

314 [-.5 * (7. / 9.)**.5, .5 * (63. / 143)**.5, -.5 * 18. / 143.**.5]]

315

316

317class VNonLocal:

318 """Wrapper class for one nonlocal term of an HGH potential."""

319 def __init__(self, l, r0, h_n):

320 self.l = l

321 self.r0 = r0

322 h_n = np.array(h_n)

323 nn = len(h_n)

324 self.nn = nn

325 self.h_n = h_n

326

327 def expand_hamiltonian_diagonal(self):

328 """Construct full atomic Hamiltonian from diagonal elements."""

329 nn = self.nn

330 h_n = self.h_n

331 h_nn = np.zeros((nn, nn))

332 for n, h in enumerate(h_n):

333 h_nn[n, n] = h

334 if self.l > 2:

335 # print 'Warning: no diagonal elements for l=%d' % l

336 # Some elements have projectors corresponding to l=3, but

337 # the HGH article only specifies how to calculate the

338 # diagonal elements of the atomic hamiltonian for l = 0, 1, 2 !

339 return

340 coefs = hcoefs_l[self.l]

341 if nn > 2:

342 h_nn[0, 2] = h_nn[2, 0] = coefs[1] * h_n[2]

343 h_nn[1, 2] = h_nn[2, 1] = coefs[2] * h_n[2]

344 if nn > 1:

345 h_nn[0, 1] = h_nn[1, 0] = coefs[0] * h_n[1]

346 return h_nn

347

348 def copy(self):

349 return VNonLocal(self.l, self.r0, self.h_n.copy())

350

351 def serialize(self): # no spin-orbit part

352 return ' '.join([' ', '%-10s' % self.r0] +

353 ['%10f' % h for h in self.h_n])

354

355

356class HGHParameterSet:

357 """Wrapper class for HGH-specific data corresponding to one element."""

358 def __init__(self, symbol, Z, Nv, rloc, c_n, v_l):

359 self.symbol = symbol # Identifier, e.g. 'Na', 'Na.sc', ...

360 self.Z = Z # Actual atomic number

361 self.Nv = Nv # Valence electron count

362 self.rloc = rloc # Characteristic radius of local part

363 self.c_n = np.array(c_n) # Polynomial coefficients for local part

364 self.v_l = list(v_l) # Non-local parts

365

366 Z, nlfe_j = configurations[self.symbol.split('.')[0]]

367 self.configuration = nlfe_j

368

369 def __str__(self):

370 strings = ['HGH setup for %s\n' % self.symbol,

371 ' Valence Z=%d, rloc=%.05f\n' % (self.Nv, self.rloc)]

372

373 if len(self.c_n) > 0:

374 coef_string = ', '.join(['%.05f' % c for c in self.c_n])

375 else:

376 coef_string = 'zeros'

377 strings.append(' Local part coeffs: %s\n' % coef_string)

378 strings.append(' Projectors:\n')

379 if not self.v_l:

380 strings.append(' None\n')

381 for v in self.v_l:

382 strings.append(' l=%d, rc=%.05f\n' % (v.l, v.r0))

383 strings.append(' Diagonal coefficients of nonlocal parts:')

384 if not self.v_l:

385 strings.append('\n None\n')

386 for v in self.v_l:

387 strings.append('\n')

388 strings.append(' l=%d: ' % v.l +

389 ', '.join(['%8.05f' % h for h in v.h_n]))

390 return ''.join(strings)

391

392 def copy(self):

393 other = HGHParameterSet(self.symbol, self.Z, self.Nv, self.rloc,

394 self.c_n, self.v_l)

395 return other

396

397 def print_info(self, txt):

398 txt(str(self))

399 txt()

400

401 def nl_iter(self):

402 for n in range(4):

403 for l, v in enumerate(self.v_l):

404 if n < v.nn:

405 yield n, l

406

407 def get_occupation_numbers(self):

408 nlfe_j = list(self.configuration)

409 nlfe_j.reverse()

410 f_ln = [[], [], []] # [[s], [p], [d]]

411 # f states will be ignored as the atomic Hamiltonians

412 # of those are, carelessly, not defined in the article.

413 lmax = len(self.v_l) - 1

414 Nv = 0

415 # Right. We need to find the occupation numbers of each state and

416 # put them into a nice list of lists f_ln.

417 #

418 # We loop over states starting with the least bound one

419 # (i.e. reversed nlfe_j), adding the occupation numbers of each state

420 # as appropriate. Once we have the right number of electrons, we

421 # end the loop.

422 #

423 # Some states in the standard configuration might

424 # be f-type; these should be skipped (unless the HGH setup actually

425 # has a valence f-state; however as noted above, some of the

426 # parameters are undefined in that case so are ignored anyway). More

427 # generally if for some state l > lmax,

428 # we can skip that state.

429 for n, l, f, e in nlfe_j:

430 if l > lmax:

431 continue

432 Nv += f

433 f_n = f_ln[l]

434 assert f_n == [] or self.symbol.endswith('.sc')

435 f_n.append(f)

436 if Nv >= self.Nv:

437 break

438 assert Nv == self.Nv

439 return f_ln

440

441 def zeropad(self):

442 """Return a new HGHParameterSet with all arrays zero padded so they

443 have the same (max) length for all such HGH setups. Makes

444 plotting multiple HGH setups easier because they have compatible

445 arrays."""

446 c_n = np.zeros(4)

447 for n, c in enumerate(self.c_n):

448 c_n[n] = c

449 v_l = []

450 for l, v in enumerate(self.v_l):

451 h_n = np.zeros(3)

452 h_n[:len(v.h_n)] = list(v.h_n)

453 v2 = VNonLocal(l, v.r0, h_n)

454 v_l.append(v2)

455 for l in range(len(self.v_l), 3):

456 v_l.append(VNonLocal(l, 0.5, np.zeros(3)))

457 copy = HGHParameterSet(self.symbol, self.Z, self.Nv, self.rloc, c_n,

458 v_l)

459 return copy

460

461 def serialize(self):

462 string1 = '%-5s %-12s %10s ' % (self.symbol, self.Z, self.rloc)

463 string2 = ' '.join(['%.10s' % c for c in self.c_n])

464 nonlocal_strings = [v.serialize() for v in self.v_l]

465 return '\n'.join([string1 + string2] + nonlocal_strings)

466

467

468def parse_local_part(string):

469 """Create HGHParameterSet object with local part initialized."""

470 tokens = iter(string.split())

471 symbol = next(tokens)

472 actual_chemical_symbol = symbol.split('.')[0]

473 Z = atomic_numbers[actual_chemical_symbol]

474 Nv = int(next(tokens))

475 rloc = float(next(tokens))

476 c_n = [float(token) for token in tokens]

477 return symbol, Z, Nv, rloc, c_n

478

479

480class HGHBogusNumbersError(ValueError):

481 """Error which is raised when the HGH parameters contain f-type

482 or higher projectors. The HGH article only defines atomic Hamiltonian

483 matrices up to l=2, so these are meaningless."""

484 pass

485

486

487def parse_hgh_setup(lines):

488 """Initialize HGHParameterSet object from text representation."""

489 lines = iter(lines)

490 symbol, Z, Nv, rloc, c_n = parse_local_part(next(lines))

491

492 def pair_up_nonlocal_lines(lines):

493 try:

494 yield (next(lines), '')

495 while True:

496 yield (next(lines), next(lines))

497 except StopIteration:

498 return

499

500 v_l = []

501 for l, (non_local, spinorbit) in enumerate(pair_up_nonlocal_lines(lines)):

502 # we discard the spinorbit 'k_n' data so far

503 nltokens = non_local.split()

504 r0 = float(nltokens[0])

505 h_n = [float(token) for token in nltokens[1:]]

506

507 # if h_n[-1] == 0.0: # Only spin-orbit contributes. Discard.

508 # h_n.pop()

509 # Actually the above causes trouble. Probably it messes up state

510 # ordering or something else that shouldn't have any effect.

511

512 vnl = VNonLocal(l, r0, h_n)

513 v_l.append(vnl)

514 if l > 2:

515 raise HGHBogusNumbersError

516

517 hgh = HGHParameterSet(symbol, Z, Nv, rloc, c_n, v_l)

518 return hgh

519

520

521def str2hgh(string):

522 return parse_hgh_setup(string.splitlines())

523

524

525def hgh2str(hgh):

526 return hgh.serialize()

527

528

529def parse_setups(lines):

530 """Read HGH data from file."""

531 setups = {}

532 entry_lines = [i for i in range(len(lines))

533 if lines[i][0].isalpha()]

534 lines_by_element = [lines[entry_lines[i]:entry_lines[i + 1]]

535 for i in range(len(entry_lines) - 1)]

536 lines_by_element.append(lines[entry_lines[-1]:])

537

538 for elines in lines_by_element:

539 try:

540 hgh = parse_hgh_setup(elines)

541 except HGHBogusNumbersError:

542 continue

543 assert hgh.symbol not in setups

544 setups[hgh.symbol] = hgh

545 return setups

546

547

548def plot(symbol, extension=None):

549 import matplotlib.pyplot as plt

550 try:

551 s = HGHSetupData(symbol)

552 except IndexError:

553 print('Nooooo')

554 return

555 s.plot()

556 if extension is not None:

557 plt.savefig(f'hgh.{symbol}.{extension}')

558

559

560def plot_many(*symbols):

561 import matplotlib.pyplot as plt

562 if not symbols:

563 symbols = setups.keys() + [key + '.sc' for key in sc_setups.keys()]

564 for symbol in symbols:

565 plt.figure(1)

566 plot(symbol, extension='png')

567 plt.clf()

568

569

570def parse_default_setups():

571 from gpaw.hgh_parameters import parameters

572 lines = parameters.splitlines()

573 setups0 = parse_setups(lines)

574 for key, value in setups0.items():

575 if key.endswith('.sc'):

576 sym, sc = key.split('.')

577 sc_setups[sym] = value

578 else:

579 setups[key] = value

580

581

582parse_default_setups()