Coverage for gpaw/atom/radialgd.py: 79%

1import numbers

2from abc import ABC, abstractmethod

3from math import factorial as fac

4from math import pi

5from typing import Tuple

7import numpy as np

8from gpaw.sphere.integrate import integrate_radial_grid

9from gpaw.spline import Spline

10from gpaw.typing import Array1D

11from gpaw.utilities import divrl, hartree

12from scipy.interpolate import make_interp_spline, splder

15def radial_grid_descriptor(eq: str, **kwargs) -> 'RadialGridDescriptor':

16 if eq == 'r=d*i':

17 assert int(kwargs['istart']) == 0

18 return EquidistantRadialGridDescriptor(float(kwargs['d']),

19 int(kwargs['iend']) + 1)

20 if eq == 'r=a*i/(n-i)':

21 beta = float(kwargs['a'])

22 ng = int(kwargs['n'])

23 return AERadialGridDescriptor(beta / ng, 1.0 / ng, ng)

24 if eq == 'r=a*i/(1-b*i)':

25 a = float(kwargs['a'])

26 b = float(kwargs['b'])

27 N = int(kwargs['n'])

28 return AERadialGridDescriptor(a, b, N)

29 if eq == 'r=a*(exp(d*i)-1)':

30 a = float(kwargs['a'])

31 d = float(kwargs['d'])

32 N = int(kwargs['iend']) + 1

33 return AbinitRadialGridDescriptor(a, d, N)

34 raise ValueError('Unknown grid: ' + eq)

37def fsbt(l, f_g, r_g, G_k):

38 """Fast spherical Bessel transform.

40 Returns::

42 oo

43 / 2

44 |r dr j (Gr) f(r),

45 / l

46 0

48 using l+1 fft's."""

50 N = (len(G_k) - 1) * 2

51 f_k = 0.0

52 F_g = f_g * r_g

53 for n in range(l + 1):

54 f_k += (r_g[1] * (1j)**(l + 1 - n) *

55 fac(l + n) / fac(l - n) / fac(n) / 2**n *

56 np.fft.rfft(F_g, N)).real * G_k**(l - n)

57 F_g[1:] /= r_g[1:]

59 f_k[1:] /= G_k[1:]**(l + 1)

60 if l == 0:

61 f_k[0] = np.dot(r_g, f_g * r_g) * r_g[1]

62 return f_k

65class RadialGridDescriptor(ABC):

66 ndim = 1 # dimension of ndarrays

68 def __init__(self, r_g: np.ndarray, dr_g, default_spline_points=25):

69 """Grid descriptor for radial grid."""

70 self.r_g = r_g

71 self.dr_g = dr_g

72 self.N = len(r_g)

73 self.dv_g = 4 * pi * r_g**2 * dr_g

74 self.default_spline_points = default_spline_points

76 def __len__(self):

77 return self.N

79 @classmethod

80 def new(cls, name, *args, **kwargs):

81 if name == 'equidistant':

82 return EquidistantRadialGridDescriptor(*args, **kwargs)

83 raise ValueError(f'Unknown grid-type: {name}')

85 def zeros(self, x=()):

86 a_xg = self.empty(x)

87 a_xg[:] = 0

88 return a_xg

90 def empty(self, x=()):

91 if isinstance(x, int):

92 x = (x,)

93 return np.zeros(x + (self.N,))

95 def integrate(self, a_xg, n=0):

96 assert n >= -2

97 return np.dot(a_xg[..., 1:],

98 (self.r_g**(2 + n) * self.dr_g)[1:]) * (4 * pi)

100 def integrate_trapz(self, a_xg, rcut=None):

101 return integrate_radial_grid(a_xg, self.r_g, rcut=rcut)

102

103 def yukawa(self, n_g, l=0, gamma=1e-6):

104 r"""Calculates the radial grid yukawa integral.

105

106 The the integral kernel for the Yukawa interaction:

107

108 \ _ _

109 exp(- /\ | r - r' |)

110 ----------------------

111 _ _

112 | r - r' |

113

114 is defined as

115

116 __ __ \ r \ r * ^ ^

117 \ 4 || I_(l+0.5)(/\ <) K_(l+0.5) (/\ >) Y (r) Y(r')

118 ) -------------------------- lm lm

119 / __ (rr')^0.5

120 lm

121

122 where I and K are the modified Bessel functions of the first

123 and second kind (K is also known as Macdonald function).

124 r = min (r, r') r = max(r, r')

125 < >

126 We now calculate the integral:

127

128

129 ^ / _ ^

130 v (r) Y (r) = |dr' n(r') Y (r')

131 l lm / l lm

132

133 with the Yukawa kernel mentioned above.

134

135 And the output array is 'vr' as it is

136 within the Hartree / radial Poisson solver.

137 """

138

139 from scipy.special import iv, kv

140 vr_g = self.zeros()

141 nrdr_g = n_g * self.r_g**1.5 * self.dr_g

142 p = 0

143 q = 0

144 k_rgamma = kv(l + 0.5, self.r_g * gamma) # K(>)

145 i_rgamma = iv(l + 0.5, self.r_g * gamma) # I(<)

146 k_rgamma[0] = kv(l + 0.5, self.r_g[1] * gamma * 1e-5)

147 # We have two integrals: one for r< and one for r>

148 # This loop-technique helps calculate them in once

149 for g_ind in range(len(nrdr_g) - 1, -1, -1):

150 dp = k_rgamma[g_ind] * nrdr_g[g_ind] # r' is r>

151 dq = i_rgamma[g_ind] * nrdr_g[g_ind] # r' is r<

152 vr_g[g_ind] = (p + 0.5 * dp) * i_rgamma[g_ind] - \

153 (q + 0.5 * dq) * k_rgamma[g_ind]

154 p += dp

155 q += dq

156 vr_g[:] += q * k_rgamma[:]

157 vr_g *= 4 * pi

158 vr_g[:] *= self.r_g[:]**0.5

159 return vr_g

160

161 def derivative_spline(self, n_g, dndr_g=None, *, degree=5):

162 """Spline-based derivative of radial function."""

163 if dndr_g is None:

164 dndr_g = self.empty()

165

166 # Boundary conditions

167 assert degree in [3, 5]

168 if degree == 5:

169 bc_type = ([(2, 0.0), (4, 0.0)], [(2, 0.0), (4, 0.0)])

170 elif degree == 3:

171 bc_type = ([(2, 0.0)], [(2, 0.0)])

172

173 s = make_interp_spline(self.r_g, n_g, k=degree, bc_type=bc_type)

174 s = splder(s)

175 dndr_g[:] = s(self.r_g)

176 return dndr_g

177

178 def derivative(self, n_g, dndr_g=None):

179 """Finite-difference derivative of radial function."""

180 if dndr_g is None:

181 dndr_g = self.empty()

182 dndr_g[0] = n_g[1] - n_g[0]

183 dndr_g[1:-1] = 0.5 * (n_g[2:] - n_g[:-2])

184 dndr_g[-1] = n_g[-1] - n_g[-2]

185 dndr_g /= self.dr_g

186 return dndr_g

187

188 def derivative2(self, a_g, b_g):

189 """Finite-difference derivative of radial function.

190

191 For an infinitely dense grid, this method would be identical

192 to the `derivative` method."""

193

194 c_g = a_g / self.dr_g

195 b_g[0] = 0.5 * c_g[1] + c_g[0]

196 b_g[1] = 0.5 * c_g[2] - c_g[0]

197 b_g[1:-1] = 0.5 * (c_g[2:] - c_g[:-2])

198 b_g[-2] = c_g[-1] - 0.5 * c_g[-3]

199 b_g[-1] = -c_g[-1] - 0.5 * c_g[-2]

200

201 def laplace(self, n_g, d2ndr2_g=None):

202 """Laplace of radial function."""

203 if d2ndr2_g is None:

204 d2ndr2_g = self.empty()

205 dndg_g = 0.5 * (n_g[2:] - n_g[:-2])

206 d2ndg2_g = n_g[2:] - 2 * n_g[1:-1] + n_g[:-2]

207 d2ndr2_g[1:-1] = (d2ndg2_g / self.dr_g[1:-1]**2 +

208 dndg_g * (self.d2gdr2()[1:-1] +

209 2 / self.r_g[1:-1] / self.dr_g[1:-1]))

210 d2ndr2_g[0] = d2ndr2_g[1]

211 d2ndr2_g[-1] = d2ndr2_g[-2]

212 return d2ndr2_g

213

214 def T(self, u_g, l):

215 dudg_g = 0.5 * (u_g[2:] - u_g[:-2])

216 d2udg2_g = u_g[2:] - 2 * u_g[1:-1] + u_g[:-2]

217 Tu_g = self.empty()

218 Tu_g[1:-1] = -0.5 * (d2udg2_g / self.dr_g[1:-1]**2 +

219 dudg_g * self.d2gdr2()[1:-1])

220 Tu_g[-1] = Tu_g[-2]

221 Tu_g[1:] += 0.5 * l * (l + 1) * u_g[1:] / self.r_g[1:]**2

222 Tu_g[0] = Tu_g[1]

223 return Tu_g

224

225 def interpolate(self, f_g, r_x):

226 from scipy.interpolate import InterpolatedUnivariateSpline

227 return InterpolatedUnivariateSpline(self.r_g, f_g)(r_x)

228

229 def fft(self, fr_g, l=0, N=None):

230 """Fourier transform.

231

232 Returns G and f(G) arrays::

233

234 _ _

235 l ^ / _ ^ iG.r

236 f(G)i Y (G) = |dr f(r)Y (r) e .

237 lm / lm

238 """

239

240 if N is None:

241 N = 2**13

242

243 assert N % 2 == 0

244

245 r_x = np.linspace(0, self.r_g[-1], N)

246 f_x = self.interpolate(fr_g, r_x)

247 f_x[1:] /= r_x[1:]

248 f_x[0] = f_x[1]

249 G_k = np.linspace(0, pi / r_x[1], N // 2 + 1)

250 f_k = 4 * pi * fsbt(l, f_x, r_x, G_k)

251 return G_k, f_k

252

253 def filter(self, f_g, rcut, Gcut, l=0):

254 Rcut = 100.0

255 N = 1024 * 8

256 r_x = np.linspace(0, Rcut, N, endpoint=False)

257 h = Rcut / N

258

259 alpha = 4.0 / rcut**2

260 mcut = np.exp(-alpha * rcut**2)

261 r2_x = r_x**2

262 m_x = np.exp(-alpha * r2_x)

263 for n in range(2):

264 m_x -= (alpha * (rcut**2 - r2_x))**n * (mcut / fac(n))

265 xcut = int(np.ceil(rcut / r_x[1]))

266 m_x[xcut:] = 0.0

267

268 G_k = np.linspace(0, pi / h, N // 2 + 1)

269

270 # Zeropad the function to same length as coordinates:

271 fpad_g = np.zeros(len(self.r_g))

272 fpad_g[:len(f_g)] = f_g

273 f_g = fpad_g

274

275 from scipy.interpolate import InterpolatedUnivariateSpline

276 if l < 2:

277 f_x = InterpolatedUnivariateSpline(self.r_g, f_g)(r_x)

278 else:

279 a_g = f_g.copy()

280 a_g[1:] /= self.r_g[1:]**(l - 1)

281 f_x = InterpolatedUnivariateSpline(

282 self.r_g, a_g)(r_x) * r_x**(l - 1)

283

284 f_x[:xcut] /= m_x[:xcut]

285 f_k = fsbt(l, f_x, r_x, G_k)

286 kcut = int(Gcut / G_k[1])

287 f_k[kcut:] = 0.0

288 ff_x = fsbt(l, f_k, G_k, r_x[:N // 2 + 1]) / pi * 2

289 ff_x *= m_x[:N // 2 + 1]

290

291 if l < 2:

292 f_g = InterpolatedUnivariateSpline(

293 r_x[:xcut + 1], ff_x[:xcut + 1])(self.r_g)

294 else:

295 ff_x[1:xcut + 1] /= r_x[1:xcut + 1]**(l - 1)

296 f_g = InterpolatedUnivariateSpline(

297 r_x[:xcut + 1], ff_x[:xcut + 1])(self.r_g) * self.r_g**(l - 1)

298 f_g[self.ceil(rcut):] = 0.0

299

300 return f_g

301

302 def poisson(self, n_g, l=0):

303 vr_g = self.zeros()

304 nrdr_g = n_g * self.r_g * self.dr_g

305 hartree(l, nrdr_g, self.r_g, vr_g)

306 return vr_g

307

308 def pseudize(self, a_g, gc, l=0, points=4):

309 """Construct smooth continuation of a_g for g<gc.

310

311 Returns (b_g, c_p[P-1]) such that b_g=a_g for g >= gc and::

312

313 P-1 2(P-1-p)+l

314 b = Sum c_p r

315 g p=0 g

316

317 for g < gc+P.

318 """

319 assert isinstance(gc, numbers.Integral) and gc > 10, gc

320

321 r_g = self.r_g

322 i = np.arange(gc, gc + points)

323 r_i = r_g[i]

324 c_p = np.polyfit(r_i**2, a_g[i] / r_i**l, points - 1)

325 b_g = a_g.copy()

326 b_g[:gc] = np.polyval(c_p, r_g[:gc]**2) * r_g[:gc]**l

327 return b_g, c_p[-1]

328

329 def cut(self, a_g, rcut):

330 gcut = self.floor(rcut)

331 r0 = 0.7 * rcut

332 x_g = np.clip((self.r_g - r0) / (rcut - r0), 0, 1)

333 f_g = x_g**2 * (3 - 2 * x_g)

334 shift = (4 * a_g[gcut] - a_g[gcut - 1]) / 3

335 a_g -= f_g * shift

336 a_g[gcut + 1:] = 0

337

338 def pseudize_normalized(self, a_g, gc, l=0, points=3):

339 """Construct normalized smooth continuation of a_g for g<gc.

340

341 Same as pseudize() with also this constraint::

342

343 / _ 2 / _ 2

344 | dr b(r) = | dr a(r)

345 / /

346 """

347

348 b_g = self.pseudize(a_g, gc, l, points)[0]

349 c_x = np.empty(points + 1)

350 gc0 = gc // 2

351 x0 = b_g[gc0]

352 r_g = self.r_g

353 i = [gc0] + list(range(gc, gc + points))

354 r_i = r_g[i]

355 norm = self.integrate(a_g**2)

356

357 def f(x):

358 b_g[gc0] = x[0]

359 c_x[:] = np.polyfit(r_i**2, b_g[i] / r_i**l, points)

360 b_g[:gc] = np.polyval(c_x, r_g[:gc]**2) * r_g[:gc]**l

361 return self.integrate(b_g**2) - norm

362

363 from scipy.optimize import fsolve

364 fsolve(f, x0)

365 return b_g, c_x[-1]

366

367 def pseudize_smooth(self,

368 a_g: Array1D,

369 gc: int,

370 l: int = 0,

371 points: int = 3,

372 Gcut: float = 6.0) -> Tuple[Array1D, float]:

373 """Minimize Fourier components above Gcut."""

374 b_g, _ = self.pseudize(a_g, gc, l, points)

375 c_x = np.empty(points + 1)

376 gc0 = gc // 2

377 x0 = b_g[gc0]

378 i = [gc0] + list(range(gc, gc + points))

379 r_g = self.r_g

380 r_i = r_g[i]

381

382 def f(x):

383 b_g[gc0] = x[0]

384 c_x[:] = np.polyfit(r_i**2, b_g[i] / r_i**l, points)

385 b_g[:gc] = np.polyval(c_x, r_g[:gc]**2) * r_g[:gc]**l

386 g_k, b_k = self.fft(b_g * r_g, l)

387 kc = int(Gcut / g_k[1])

388 f_k = g_k[kc:] * b_k[kc:]

389 return f_k @ f_k

390

391 from scipy.optimize import minimize

392 minimize(f, x0)

393

394 return b_g, c_x[-1]

395

396 def jpseudize(self, a_g, gc, l=0, points=4):

397 """Construct spherical Bessel function continuation of a_g for g<gc.

398

399 Returns (b_g, b(0)/r^l) such that b_g=a_g for g >= gc and::

400

401 P-2

402 b = Sum c_p j (q r )

403 g p=0 l p g

404

405 for g < gc+P, where.

406 """

407

408 from scipy.optimize import brentq

409 from scipy.special import sph_jn

410

411 if a_g[gc] == 0:

412 return self.zeros(), 0.0

413

414 assert isinstance(gc, int) and gc > 10

415

416 zeros_l = [[1, 2, 3, 4, 5, 6],

417 [1.430, 2.459, 3.471, 4.477, 5.482, 6.484],

418 [1.835, 2.895, 3.923, 4.938, 5.949, 6.956],

419 [2.224, 3.316, 4.360, 5.387, 6.405, 7.418]]

420

421 # Logarithmic derivative:

422 ld = np.dot([-1 / 60, 3 / 20, -3 / 4, 0, 3 / 4, -3 / 20, 1 / 60],

423 a_g[gc - 3:gc + 4]) / a_g[gc] / self.dr_g[gc]

424

425 rc = self.r_g[gc]

426

427 def f(q):

428 j, dj = (y[-1] for y in sph_jn(l, q * rc))

429 return dj * q - j * ld

430

431 j_pg = self.empty(points - 1)

432 q_p = np.empty(points - 1)

433

434 zeros = zeros_l[l]

435 if rc * ld > l:

436 z1 = zeros[0]

437 zeros = zeros[1:]

438 else:

439 z1 = 0

440 x = 0.01

441 for p, z2 in enumerate(zeros[:points - 1]):

442 q = brentq(f, z1 * pi / rc + x, z2 * pi / rc - x)

443 j_pg[p] = [sph_jn(l, q * r)[0][-1] for r in self.r_g]

444 q_p[p] = q

445 z1 = z2

446

447 C_dg = [[0, 0, 0, 1, 0, 0, 0],

448 [1 / 90, -3 / 20, 3 / 2, -49 / 18, 3 / 2, -3 / 20, 1 / 90],

449 [1 / 8, -1, 13 / 8, 0, -13 / 8, 1, -1 / 8],

450 [-1 / 6, 2, -13 / 2, 28 / 3, -13 / 2, 2, -1 / 6]][:points - 1]

451 c_p = np.linalg.solve(np.dot(C_dg, j_pg[:, gc - 3:gc + 4].T),

452 np.dot(C_dg, a_g[gc - 3:gc + 4]))

453 b_g = a_g.copy()

454 b_g[:gc + 2] = np.dot(c_p, j_pg[:, :gc + 2])

455 return b_g, np.dot(c_p, q_p**l) * 2**l * fac(l) / fac(2 * l + 1)

456

457 def plot(self, a_g, n=0, rc=4.0, show=False):

458 import matplotlib.pyplot as plt

459 r_g = self.r_g[:len(a_g)]

460 if n < 0:

461 r_g = r_g[1:]

462 a_g = a_g[1:]

463 plt.plot(r_g, a_g * r_g**n)

464 plt.axis(xmin=0, xmax=rc)

465 if show:

466 plt.show()

467

468 def floor(self, r):

469 return np.floor(self.r2g(r)).astype(int)

470

471 def round(self, r):

472 return np.around(self.r2g(r)).astype(int)

473

474 def ceil(self, r):

475 return np.ceil(self.r2g(r)).astype(int)

476

477 def spline(self, a_g, rcut=None, l=0, points=None,

478 backwards_compatible=True):

479 """Create spline representation of a radial function f(r).

480

481 The spline represents a rescaled version of the function, f(r) / r^l.

482 """

483 if points is None:

484 points = self.default_spline_points

485

486 if rcut is None:

487 # Calculate rcut for the input function, i.e. a value rcut which

488 # satisfies f(r) = 0 ∀ r ≥ rcut

489 g = len(a_g) - 1

490 while a_g[g] == 0.0:

491 g -= 1

492 rcut = self.r_g[g + 1]

493

494 # Rescale f(r) -> f(r) / r^l

495 N = len(a_g)

496 b_g = divrl(a_g, l, self.r_g[:N])

497

498 # Interpolate to a uniform radial grid (for the spline representation)

499 r_i = np.linspace(0, rcut, points + 1)

500 if not backwards_compatible:

501 from scipy.interpolate import CubicSpline

502 b_i = CubicSpline(self.r_g[:N], b_g, bc_type='clamped')(r_i)

503 else:

504 g_i = np.clip((self.r2g(r_i) + 0.5).astype(int), 1, N - 2)

505 r1_i = self.r_g[g_i - 1]

506 r2_i = self.r_g[g_i]

507 r3_i = self.r_g[g_i + 1]

508 x1_i = (r_i - r2_i) * (r_i - r3_i) / (r1_i - r2_i) / (r1_i - r3_i)

509 x2_i = (r_i - r1_i) * (r_i - r3_i) / (r2_i - r1_i) / (r2_i - r3_i)

510 x3_i = (r_i - r1_i) * (r_i - r2_i) / (r3_i - r1_i) / (r3_i - r2_i)

511 b1_i = b_g[g_i - 1]

512 b2_i = b_g[g_i]

513 b3_i = b_g[g_i + 1]

514 b_i = b1_i * x1_i + b2_i * x2_i + b3_i * x3_i

515

516 return Spline.from_data(l, rcut, b_i)

517

518 def get_cutoff(self, f_g):

519 g = self.N - 1

520 while f_g[g] == 0.0:

521 g -= 1

522 gcut = g + 1

523 return gcut

524

525 @abstractmethod

526 def r2g(self, r):

527 """Inverse continuous map from a real space distance (r)

528 to a floating point grid index (g).

529

530 Used by methods floor, round, and ceil, which manipulate this

531 floating point to an integer accordingly.

532 """

533

534

535class EquidistantRadialGridDescriptor(RadialGridDescriptor):

536 def __init__(self, h, N=1000, h0=0.0):

537 """Equidistant radial grid descriptor.

538

539 The radial grid is r(g) = h0 + g*h, g = 0, 1, ..., N - 1."""

540

541 RadialGridDescriptor.__init__(self,

542 h * np.arange(N) + h0,

543 h + np.zeros(N))

544

545 def r2g(self, r):

546 return int(np.ceil((r - self.r_g[0]) / (self.r_g[1] - self.r_g[0])))

547

548 def xml(self, id='grid1'):

549 assert self.r_g[0] == 0.0

550 return ('<radial_grid eq="r=d*i" d="{}" '

551 'istart="0" iend="{}" id="{}"/>\n'

552 .format(self.r_g[1], len(self.r_g) - 1, id))

553

554 def spline(self, a_g, rcut=None, l=0, points=None):

555 b_g = a_g.copy()

556 if l > 0:

557 b_g = divrl(b_g, l, self.r_g[:len(a_g)])

558 return Spline.from_data(l, self.r_g[len(a_g) - 1], b_g)

559

560

561class AERadialGridDescriptor(RadialGridDescriptor):

562 def __init__(self, a, b, N=1000, default_spline_points=25):

563 """Radial grid descriptor for all-electron calculation.

564

565 The radial grid is::

566

567 a g

568 r(g) = -------, g = 0, 1, ..., N - 1

569 1 - b g

570 """

571

572 self.a = a

573 self.b = b

574 g = np.arange(N)

575 r_g = self.a * g / (1 - self.b * g)

576 dr_g = (self.b * r_g + self.a)**2 / self.a

577 RadialGridDescriptor.__init__(self, r_g, dr_g, default_spline_points)

578 self._d2gdr2 = -2 * self.a * self.b / (self.b * self.r_g + self.a)**3

579

580 @property

581 def beta(self):

582 return self.a * self.N

583

584 def __repr__(self):

585 return (

586 f'{type(self).__name__}'

587 f'({self.a}, {self.b}, {self.N}, {self.default_spline_points})')

588

589 def r2g(self, r):

590 # return r / (r * self.b + self.a)

591 # Hack to preserve backwards compatibility:

592 ng = 1.0 / self.b

593 beta = self.a / self.b

594 return ng * r / (beta + r)

595

596 def new(self, N):

597 return AERadialGridDescriptor(self.a, self.b, N)

598

599 def xml(self, id='grid1'):

600 if abs(self.N - 1 / self.b) < 1e-5:

601 return (('<radial_grid eq="r=a*i/(n-i)" a="%r" n="%d" ' +

602 'istart="0" iend="%d" id="%s"/>\n') %

603 (self.a * self.N, self.N, self.N - 1, id))

604 return (('<radial_grid eq="r=a*i/(1-b*i)" a="%r" b="%r" n="%d" ' +

605 'istart="0" iend="%d" id="%s"/>\n') %

606 (self.a, self.b, self.N, self.N - 1, id))

607

608 def d2gdr2(self):

609 return self._d2gdr2

610

611

612class AbinitRadialGridDescriptor(RadialGridDescriptor):

613 def __init__(self, a, d, N=1000, default_spline_points=25):

614 """Radial grid descriptor for Abinit calculations.

615

616 The radial grid is::

617

618 dg

619 r(g) = a(e - 1), g = 0, 1, ..., N - 1

620 """

621

622 self.a = a

623 self.d = d

624 g = np.arange(N)

625 r_g = a * (np.exp(d * g) - 1)

626 dr_g = (r_g + a) * d

627 RadialGridDescriptor.__init__(self, r_g, dr_g, default_spline_points)

628

629 def r2g(self, r):

630 return np.log(r / self.a + 1) / self.d

631

632 def d2gdr2(self):

633 return -1 / (self.a**2 * self.d * (self.r_g / self.a + 1)**2)

634

635 def new(self, N):

636 return AbinitRadialGridDescriptor(self.a, self.d, N)