Coverage for gpaw/atom/generator.py: 78%

1# Copyright (C) 2003 CAMP 

2# Please see the accompanying LICENSE file for further information. 

3from functools import cached_property 

4from math import pi, sqrt 

5 

6import numpy as np 

7from numpy.linalg import solve, inv 

8from scipy.linalg import eigh 

9from ase.units import Ha 

10 

11from gpaw.setup_data import SetupData 

12from gpaw.atom.configurations import configurations 

13from gpaw import __version__ as version 

14from gpaw.atom.all_electron import AllElectron, shoot 

15from gpaw.utilities import hartree 

16from gpaw.xc.hybrid import constructX, atomic_exact_exchange 

17from gpaw.atom.filter import Filter 

18 

19 

20def enumerate_j_ln(anything_ln: list[list[object]]) -> list[list[int]]: 

21 """Return j_ln to map between _j and _ln index.""" 

22 return [[0 for _ in anything_n] for anything_n in anything_ln] 

23 

24 

25class Generator(AllElectron): 

26 def __init__(self, symbol, xcname='LDA', scalarrel=True, corehole=None, 

27 configuration=None, 

28 nofiles=True, txt='-', gpernode=150, 

29 orbital_free=False, tw_coeff=1.): 

30 AllElectron.__init__(self, symbol, xcname, scalarrel, corehole, 

31 configuration, nofiles, txt, gpernode, 

32 orbital_free, tw_coeff) 

33 

34 def run(self, core='', rcut=1.0, extra=None, 

35 logderiv=False, vbar=None, exx=False, name=None, 

36 normconserving='', filter=(0.4, 1.75), rcutcomp=None, 

37 write_xml=True, 

38 empty_states='', yukawa_gamma=0.0): 

39 

40 self.name = name 

41 

42 self.core = core 

43 if isinstance(rcut, float): 

44 rcut_l = [rcut] 

45 else: 

46 rcut_l = rcut 

47 rcutmax = max(rcut_l) 

48 rcutmin = min(rcut_l) 

49 self.rcut_l = rcut_l 

50 

51 if rcutcomp is None: 

52 rcutcomp = rcutmin 

53 self.rcutcomp = rcutcomp 

54 

55 hfilter, xfilter = filter 

56 

57 Z = self.Z 

58 

59 n_j = self.n_j 

60 l_j = self.l_j 

61 f_j = self.f_j 

62 e_j = self.e_j 

63 

64 if vbar is None: 

65 vbar = ('poly', rcutmin * 0.9) 

66 vbar_type, rcutvbar = vbar 

67 

68 normconserving_l = [x in normconserving for x in 'spdf'] 

69 

70 # Parse core string: 

71 j = 0 

72 if core.startswith('['): 

73 a, core = core.split(']') 

74 core_symbol = a[1:] 

75 j = len(configurations[core_symbol][1]) 

76 

77 while core != '': 

78 assert n_j[j] == int(core[0]) 

79 assert l_j[j] == 'spdf'.find(core[1]) 

80 if j != self.jcorehole: 

81 assert f_j[j] == 2 * (2 * l_j[j] + 1) 

82 j += 1 

83 core = core[2:] 

84 

85 njcore = j 

86 self.njcore = njcore 

87 

88 while empty_states != '': 

89 n = int(empty_states[0]) 

90 l = 'spdf'.find(empty_states[1]) 

91 assert n == 1 + l + l_j.count(l) 

92 n_j.append(n) 

93 l_j.append(l) 

94 f_j.append(0.0) 

95 e_j.append(-0.01) 

96 empty_states = empty_states[2:] 

97 

98 if 2 in l_j[njcore:]: 

99 # We have a bound valence d-state. Add bound s- and 

100 # p-states if not already there: 

101 for l in [0, 1]: 

102 if l not in l_j[njcore:]: 

103 n_j.append(1 + l + l_j.count(l)) 

104 l_j.append(l) 

105 f_j.append(0.0) 

106 e_j.append(-0.01) 

107 

108 if l_j[njcore:] == [0] and Z > 2: 

109 # We have only a bound valence s-state and we are not 

110 # hydrogen and not helium. Add bound p-state: 

111 n_j.append(n_j[njcore]) 

112 l_j.append(1) 

113 f_j.append(0.0) 

114 e_j.append(-0.01) 

115 

116 nj = len(n_j) 

117 

118 self.Nv = sum(f_j[njcore:]) 

119 self.Nc = sum(f_j[:njcore]) 

120 

121 lmaxocc = max(l_j[njcore:]) 

122 lmax = max(l_j[njcore:]) 

123 

124 # Parameters for orbital_free 

125 if self.orbital_free: 

126 self.n_j = [1] 

127 self.l_j = [0] 

128 self.f_j = [self.Z] 

129 self.e_j = [self.e_j[0]] 

130 

131 n_j = self.n_j 

132 l_j = self.l_j 

133 f_j = self.f_j 

134 e_j = self.e_j 

135 nj = len(n_j) 

136 lmax = 0 

137 lmaxocc = 0 

138 

139 # Do all-electron calculation: 

140 AllElectron.run(self) 

141 

142 # Highest occupied atomic orbital: 

143 self.emax = max(e_j) 

144 

145 N = self.N 

146 r = self.r 

147 dr = self.dr 

148 d2gdr2 = self.d2gdr2 

149 beta = self.beta 

150 

151 dv = r**2 * dr 

152 

153 t = self.text 

154 t() 

155 t('Generating PAW setup') 

156 if core != '': 

157 t('Frozen core:', core) 

158 

159 # So far - no ghost-states: 

160 self.ghost = False 

161 

162 # Calculate the kinetic energy of the core states: 

163 Ekincore = 0.0 

164 j = 0 

165 for f, e, u in zip(f_j[:njcore], e_j[:njcore], self.u_j[:njcore]): 

166 u = np.where(abs(u) < 1e-160, 0, u) # XXX Numeric! 

167 k = e - np.sum((u**2 * self.vr * dr)[1:] / r[1:]) 

168 Ekincore += f * k 

169 if j == self.jcorehole: 

170 self.Ekincorehole = k 

171 j += 1 

172 

173 # Calculate core density: 

174 if njcore == 0: 

175 nc = np.zeros(N) 

176 else: 

177 uc_j = self.u_j[:njcore] 

178 uc_j = np.where(abs(uc_j) < 1e-160, 0, uc_j) # XXX Numeric! 

179 nc = np.dot(f_j[:njcore], uc_j**2) / (4 * pi) 

180 nc[1:] /= r[1:]**2 

181 nc[0] = nc[1] 

182 

183 self.nc = nc 

184 

185 # Calculate core kinetic energy density 

186 if njcore == 0: 

187 tauc = np.zeros(N) 

188 else: 

189 tauc = self.radial_kinetic_energy_density(f_j[:njcore], 

190 l_j[:njcore], 

191 self.u_j[:njcore]) 

192 t('Kinetic energy of the core from tauc =', 

193 np.dot(tauc * r * r, dr) * 4 * pi) 

194 

195 # Order valence states with respect to angular momentum 

196 # quantum number: 

197 self.n_ln = n_ln = [] 

198 self.f_ln = f_ln = [] 

199 self.e_ln = e_ln = [] 

200 for l in range(lmax + 1): 

201 n_n = [] 

202 f_n = [] 

203 e_n = [] 

204 for j in range(njcore, nj): 

205 if l_j[j] == l: 

206 n_n.append(n_j[j]) 

207 f_n.append(f_j[j]) 

208 e_n.append(e_j[j]) 

209 n_ln.append(n_n) 

210 f_ln.append(f_n) 

211 e_ln.append(e_n) 

212 

213 # Add extra projectors: 

214 if extra is not None: 

215 if len(extra) == 0: 

216 lmaxextra = 0 

217 else: 

218 lmaxextra = max(extra.keys()) 

219 if lmaxextra > lmax: 

220 for l in range(lmax, lmaxextra): 

221 n_ln.append([]) 

222 f_ln.append([]) 

223 e_ln.append([]) 

224 lmax = lmaxextra 

225 for l in extra: 

226 nn = -1 

227 for e in extra[l]: 

228 n_ln[l].append(nn) 

229 f_ln[l].append(0.0) 

230 e_ln[l].append(e) 

231 nn -= 1 

232 else: 

233 # Automatic: 

234 

235 # Make sure we have two projectors for each occupied channel: 

236 for l in range(lmaxocc + 1): 

237 if len(n_ln[l]) < 2 and not normconserving_l[l]: 

238 # Only one - add one more: 

239 assert len(n_ln[l]) == 1 

240 n_ln[l].append(-1) 

241 f_ln[l].append(0.0) 

242 e_ln[l].append(1.0 + e_ln[l][0]) 

243 

244 if lmaxocc < 2 and lmaxocc == lmax: 

245 # Add extra projector for l = lmax + 1: 

246 n_ln.append([-1]) 

247 f_ln.append([0.0]) 

248 e_ln.append([0.0]) 

249 lmax += 1 

250 

251 self.lmax = lmax 

252 

253 rcut_l.extend([rcutmin] * (lmax + 1 - len(rcut_l))) 

254 

255 t('Cutoffs:') 

256 for rc, s in zip(rcut_l, 'spdf'): 

257 t(f'rc({s})={rc:.3f}') 

258 t('rc(vbar)=%.3f' % rcutvbar) 

259 t('rc(comp)=%.3f' % rcutcomp) 

260 t('rc(nct)=%.3f' % rcutmax) 

261 t() 

262 t('Kinetic energy of the core states: %.6f' % Ekincore) 

263 

264 # Allocate arrays: 

265 self.u_ln = u_ln = [] # phi * r 

266 self.s_ln = s_ln = [] # phi-tilde * r 

267 self.q_ln = q_ln = [] # p-tilde * r 

268 for l in range(lmax + 1): 

269 nn = len(n_ln[l]) 

270 u_ln.append(np.zeros((nn, N))) 

271 s_ln.append(np.zeros((nn, N))) 

272 q_ln.append(np.zeros((nn, N))) 

273 

274 # Fill in all-electron wave functions: 

275 for l in range(lmax + 1): 

276 # Collect all-electron wave functions: 

277 u_n = [self.u_j[j] for j in range(njcore, nj) if l_j[j] == l] 

278 for n, u in enumerate(u_n): 

279 u_ln[l][n] = u 

280 

281 # Grid-index corresponding to rcut: 

282 gcut_l = [1 + int(rc * N / (rc + beta)) for rc in rcut_l] 

283 

284 rcutfilter = xfilter * rcutmax 

285 self.rcutfilter = rcutfilter 

286 gcutfilter = 1 + int(rcutfilter * N / (rcutfilter + beta)) 

287 gcutmax = 1 + int(rcutmax * N / (rcutmax + beta)) 

288 

289 # Outward integration of unbound states stops at 3 * rcut: 

290 gmax = int(3 * rcutmax * N / (3 * rcutmax + beta)) 

291 assert gmax > gcutfilter 

292 

293 # Calculate unbound extra states: 

294 c2 = -(r / dr)**2 

295 c10 = -d2gdr2 * r**2 

296 for l, (n_n, e_n, u_n) in enumerate(zip(n_ln, e_ln, u_ln)): 

297 for n, e, u in zip(n_n, e_n, u_n): 

298 if n < 0: 

299 u[:] = 0.0 

300 shoot(u, l, self.vr, e, self.r2dvdr, r, dr, c10, c2, 

301 self.scalarrel, gmax=gmax) 

302 u *= 1.0 / u[gcut_l[l]] 

303 

304 charge = Z - self.Nv - self.Nc 

305 t('Charge: %.1f' % charge) 

306 t('Core electrons: %.1f' % self.Nc) 

307 t('Valence electrons: %.1f' % self.Nv) 

308 

309 # Construct smooth wave functions: 

310 coefs = [] 

311 for l, (u_n, s_n) in enumerate(zip(u_ln, s_ln)): 

312 nodeless = True 

313 gc = gcut_l[l] 

314 for u, s in zip(u_n, s_n): 

315 s[:] = u 

316 if normconserving_l[l]: 

317 A = np.zeros((5, 5)) 

318 A[:4, 0] = 1.0 

319 A[:4, 1] = r[gc - 2:gc + 2]**2 

320 A[:4, 2] = A[:4, 1]**2 

321 A[:4, 3] = A[:4, 1] * A[:4, 2] 

322 A[:4, 4] = A[:4, 2]**2 

323 A[4, 4] = 1.0 

324 a = u[gc - 2:gc + 3] / r[gc - 2:gc + 3]**(l + 1) 

325 a = np.log(a) 

326 

327 def f(x): 

328 a[4] = x 

329 b = solve(A, a) 

330 r1 = r[:gc] 

331 r2 = r1**2 

332 rl1 = r1**(l + 1) 

333 y = b[0] + r2 * (b[1] + r2 * (b[2] + r2 * 

334 (b[3] + r2 * b[4]))) 

335 y = np.exp(y) 

336 s[:gc] = rl1 * y 

337 return np.dot(s**2, dr) - 1 

338 x1 = 0.0 

339 x2 = 0.001 

340 f1 = f(x1) 

341 f2 = f(x2) 

342 while abs(f1) > 1e-6: 

343 x0 = (x1 / f1 - x2 / f2) / (1 / f1 - 1 / f2) 

344 f0 = f(x0) 

345 if abs(f1) < abs(f2): 

346 x2, f2 = x1, f1 

347 x1, f1 = x0, f0 

348 

349 else: 

350 A = np.ones((4, 4)) 

351 A[:, 0] = 1.0 

352 A[:, 1] = r[gc - 2:gc + 2]**2 

353 A[:, 2] = A[:, 1]**2 

354 A[:, 3] = A[:, 1] * A[:, 2] 

355 a = u[gc - 2:gc + 2] / r[gc - 2:gc + 2]**(l + 1) 

356 if 0: # l < 2 and nodeless: 

357 a = np.log(a) 

358 a = solve(A, a) 

359 r1 = r[:gc] 

360 r2 = r1**2 

361 rl1 = r1**(l + 1) 

362 y = a[0] + r2 * (a[1] + r2 * (a[2] + r2 * (a[3]))) 

363 if 0: # l < 2 and nodeless: 

364 y = np.exp(y) 

365 s[:gc] = rl1 * y 

366 

367 coefs.append(a) 

368 if nodeless: 

369 if not (s[1:gc] > 0.0).all(): 

370 raise RuntimeError( 

371 'Error: The %d%s pseudo wave has a node!' % 

372 (n_ln[l][0], 'spdf'[l])) 

373 # Only the first state for each l must be nodeless: 

374 nodeless = False 

375 

376 # Calculate pseudo core density: 

377 gcutnc = 1 + int(rcutmax * N / (rcutmax + beta)) 

378 self.nct = nct = nc.copy() 

379 A = np.ones((4, 4)) 

380 A[0] = 1.0 

381 A[1] = r[gcutnc - 2:gcutnc + 2]**2 

382 A[2] = A[1]**2 

383 A[3] = A[1] * A[2] 

384 a = nc[gcutnc - 2:gcutnc + 2] 

385 a = solve(np.transpose(A), a) 

386 r2 = r[:gcutnc]**2 

387 nct[:gcutnc] = a[0] + r2 * (a[1] + r2 * (a[2] + r2 * a[3])) 

388 t('Pseudo-core charge: %.6f' % (4 * pi * np.dot(nct, dv))) 

389 

390 # ... and the pseudo core kinetic energy density: 

391 tauct = tauc.copy() 

392 a = tauc[gcutnc - 2:gcutnc + 2] 

393 a = solve(np.transpose(A), a) 

394 tauct[:gcutnc] = a[0] + r2 * (a[1] + r2 * (a[2] + r2 * a[3])) 

395 

396 # ... and the soft valence density: 

397 nt = np.zeros(N) 

398 for f_n, s_n in zip(f_ln, s_ln): 

399 nt += np.dot(f_n, s_n**2) / (4 * pi) 

400 nt[1:] /= r[1:]**2 

401 nt[0] = nt[1] 

402 nt += nct 

403 self.nt = nt 

404 

405 # Calculate the shape function: 

406 x = r / rcutcomp 

407 gaussian = np.zeros(N) 

408 self.gamma = gamma = 10.0 

409 gaussian[:gmax] = np.exp(-gamma * x[:gmax]**2) 

410 gt = 4 * (gamma / rcutcomp**2)**1.5 / sqrt(pi) * gaussian 

411 t('Shape function alpha=%.3f' % (gamma / rcutcomp**2)) 

412 norm = np.dot(gt, dv) 

413 # print norm, norm-1 

414 assert abs(norm - 1) < 1e-2 

415 gt /= norm 

416 

417 # Calculate smooth charge density: 

418 Nt = np.dot(nt, dv) 

419 rhot = nt - (Nt + charge / (4 * pi)) * gt 

420 t('Pseudo-electron charge', 4 * pi * Nt) 

421 

422 vHt = np.zeros(N) 

423 hartree(0, rhot * r * dr, r, vHt) 

424 vHt[1:] /= r[1:] 

425 vHt[0] = vHt[1] 

426 

427 vXCt = np.zeros(N) 

428 

429 extra_xc_data = {} 

430 

431 if self.xc.type != 'GLLB': 

432 self.xc.calculate_spherical(self.rgd, 

433 nt.reshape((1, -1)), 

434 vXCt.reshape((1, -1))) 

435 else: 

436 self.xc.get_smooth_xc_potential_and_energy_1d(vXCt) 

437 

438 # Calculate extra-stuff for non-local functionals 

439 self.xc.get_extra_setup_data(extra_xc_data) 

440 

441 vt = vHt + vXCt 

442 

443 if self.orbital_free: 

444 vt /= self.tw_coeff 

445 

446 # Construct zero potential: 

447 gc = 1 + int(rcutvbar * N / (rcutvbar + beta)) 

448 if vbar_type == 'f': 

449 assert lmax == 2 

450 uf = np.zeros(N) 

451 l = 3 

452 

453 # Solve for all-electron f-state: 

454 eps = 0.0 

455 shoot(uf, l, self.vr, eps, self.r2dvdr, r, dr, c10, c2, 

456 self.scalarrel, gmax=gmax) 

457 uf *= 1.0 / uf[gc] 

458 

459 # Fit smooth pseudo f-state polynomium: 

460 A = np.ones((4, 4)) 

461 A[:, 0] = 1.0 

462 A[:, 1] = r[gc - 2:gc + 2]**2 

463 A[:, 2] = A[:, 1]**2 

464 A[:, 3] = A[:, 1] * A[:, 2] 

465 a = uf[gc - 2:gc + 2] / r[gc - 2:gc + 2]**(l + 1) 

466 a0, a1, a2, a3 = solve(A, a) 

467 r1 = r[:gc] 

468 r2 = r1**2 

469 rl1 = r1**(l + 1) 

470 y = a0 + r2 * (a1 + r2 * (a2 + r2 * a3)) 

471 sf = uf.copy() 

472 sf[:gc] = rl1 * y 

473 

474 # From 0 to gc, use analytic formula for kinetic energy operator: 

475 r4 = r2**2 

476 r6 = r4 * r2 

477 enumerator = (a0 * l * (l + 1) + 

478 a1 * (l + 2) * (l + 3) * r2 + 

479 a2 * (l + 4) * (l + 5) * r4 + 

480 a3 * (l + 6) * (l + 7) * r6) 

481 denominator = a0 + a1 * r2 + a2 * r4 + a3 * r6 

482 ekin_over_phit = - 0.5 * (enumerator / denominator - l * (l + 1)) 

483 ekin_over_phit[1:] /= r2[1:] 

484 

485 vbar = eps - vt 

486 vbar[:gc] -= ekin_over_phit 

487 vbar[0] = vbar[1] # Actually we can collect the terms into 

488 # a single fraction without poles, so as to avoid doing this, 

489 # but this is good enough 

490 

491 # From gc to gmax, use finite-difference formula for kinetic 

492 # energy operator: 

493 vbar[gc:gmax] -= self.kin(l, sf)[gc:gmax] / sf[gc:gmax] 

494 vbar[gmax:] = 0.0 

495 else: 

496 assert vbar_type == 'poly' 

497 A = np.ones((2, 2)) 

498 A[0] = 1.0 

499 A[1] = r[gc - 1:gc + 1]**2 

500 a = vt[gc - 1:gc + 1] 

501 a = solve(np.transpose(A), a) 

502 r2 = r**2 

503 vbar = a[0] + r2 * a[1] - vt 

504 vbar[gc:] = 0.0 

505 

506 vt += vbar 

507 

508 # Construct projector functions: 

509 for l, (e_n, s_n, q_n) in enumerate(zip(e_ln, s_ln, q_ln)): 

510 for e, s, q in zip(e_n, s_n, q_n): 

511 q[:] = self.kin(l, s) + (vt - e) * s 

512 q[gcutmax:] = 0.0 

513 

514 filter = Filter(r, dr, gcutfilter, hfilter).filter 

515 if self.orbital_free: 

516 vbar *= self.tw_coeff 

517 vbar = filter(vbar * r) 

518 

519 # Calculate matrix elements: 

520 self.dK_lnn = dK_lnn = [] 

521 self.dH_lnn = dH_lnn = [] 

522 self.dO_lnn = dO_lnn = [] 

523 

524 for l, (e_n, u_n, s_n, q_n) in enumerate(zip(e_ln, u_ln, 

525 s_ln, q_ln)): 

526 A_nn = np.inner(s_n, q_n * dr) 

527 # Do a LU decomposition of A: 

528 nn = len(e_n) 

529 L_nn = np.identity(nn, float) 

530 U_nn = A_nn.copy() 

531 

532 # Keep all bound states normalized 

533 if sum([n > 0 for n in n_ln[l]]) <= 1: 

534 for i in range(nn): 

535 for j in range(i + 1, nn): 

536 L_nn[j, i] = 1.0 * U_nn[j, i] / U_nn[i, i] 

537 U_nn[j, :] -= U_nn[i, :] * L_nn[j, i] 

538 

539 dO_nn = (np.inner(u_n, u_n * dr) - 

540 np.inner(s_n, s_n * dr)) 

541 

542 e_nn = np.zeros((nn, nn)) 

543 e_nn.ravel()[::nn + 1] = e_n 

544 dH_nn = np.dot(dO_nn, e_nn) - A_nn 

545 

546 q_n[:] = np.dot(inv(np.transpose(U_nn)), q_n) 

547 s_n[:] = np.dot(inv(L_nn), s_n) 

548 u_n[:] = np.dot(inv(L_nn), u_n) 

549 

550 dO_nn = np.dot(np.dot(inv(L_nn), dO_nn), 

551 inv(np.transpose(L_nn))) 

552 dH_nn = np.dot(np.dot(inv(L_nn), dH_nn), 

553 inv(np.transpose(L_nn))) 

554 

555 ku_n = [self.kin(l, u, e) for u, e in zip(u_n, e_n)] 

556 ks_n = [self.kin(l, s) for s in s_n] 

557 dK_nn = 0.5 * (np.inner(u_n, ku_n * dr) - 

558 np.inner(s_n, ks_n * dr)) 

559 dK_nn += np.transpose(dK_nn).copy() 

560 

561 dK_lnn.append(dK_nn) 

562 dO_lnn.append(dO_nn) 

563 dH_lnn.append(dH_nn) 

564 

565 for n, q in enumerate(q_n): 

566 q[:] = filter(q, l) * r**(l + 1) 

567 

568 A_nn = np.inner(s_n, q_n * dr) 

569 q_n[:] = np.dot(inv(np.transpose(A_nn)), q_n) 

570 

571 self.vt = vt 

572 self.vbar = vbar 

573 

574 t('state eigenvalue norm') 

575 t('--------------------------------') 

576 for l, (n_n, f_n, e_n) in enumerate(zip(n_ln, f_ln, e_ln)): 

577 for n in range(len(e_n)): 

578 if n_n[n] > 0: 

579 f = '(%d)' % f_n[n] 

580 t('%d%s%-4s: %12.6f %12.6f' % ( 

581 n_n[n], 'spdf'[l], f, e_n[n], 

582 np.dot(s_ln[l][n]**2, dr))) 

583 else: 

584 t('*{} : {:12.6f}'.format('spdf'[l], e_n[n])) 

585 t('--------------------------------') 

586 

587 self.logd = {} 

588 if logderiv: 

589 ni = 300 

590 self.elog = np.linspace(-5.0, 1.0, ni) 

591 delta = self.elog[1] - self.elog[0] 

592 # Calculate logarithmic derivatives: 

593 gld = gcutmax + 10 

594 self.rlog = r[gld] 

595 assert gld < gmax 

596 t('Calculating logarithmic derivatives at r=%.3f' % r[gld]) 

597 t('(skip with [Ctrl-C])') 

598 

599 try: 

600 u = np.zeros(N) 

601 for l in range(4): 

602 self.logd[l] = (np.empty(ni), np.empty(ni)) 

603 if l <= lmax: 

604 dO_nn = dO_lnn[l] 

605 dH_nn = dH_lnn[l] 

606 q_n = q_ln[l] 

607 

608 ae = self.symbol + '.ae.ld.' + 'spdf'[l] 

609 ps = self.symbol + '.ps.ld.' + 'spdf'[l] 

610 with open(ae, 'w') as fae, open(ps, 'w') as fps: 

611 ld_ae = 1000 

612 ld_ps = 1000 

613 nodes_ae = [] 

614 nodes_ps = [] 

615 for i, e in enumerate(self.elog): 

616 # All-electron logarithmic derivative: 

617 u[:] = 0.0 

618 shoot(u, l, self.vr, e, self.r2dvdr, r, dr, 

619 c10, c2, self.scalarrel, gmax=gld) 

620 dudr = 0.5 * (u[gld + 1] - u[gld - 1]) / dr[gld] 

621 ld = dudr / u[gld] - 1.0 / r[gld] 

622 if ld > ld_ae: 

623 nodes_ae.append(e) 

624 ld_ae = ld 

625 print(e, ld, file=fae) 

626 self.logd[l][0][i] = ld 

627 

628 # PAW logarithmic derivative: 

629 s = self.integrate(l, vt, e, gld) 

630 if l <= lmax: 

631 A_nn = dH_nn - e * dO_nn 

632 s_n = [self.integrate(l, vt, e, gld, q) 

633 for q in q_n] 

634 B_nn = np.inner(q_n, s_n * dr) 

635 a_n = np.dot(q_n, s * dr) 

636 

637 B_nn = np.dot(A_nn, B_nn)