Coverage for eminus/xc/lda_xc_gdsmfb.py: 100.00%

1# SPDX-FileCopyrightText: 2021 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""GDSMFB LDA exchange-correlation. 

4 

5Reference: Phys. Rev. Lett. 119, 135001. 

6""" 

7 

8import dataclasses 

9 

10import numpy as np 

11 

12 

13def lda_xc_gdsmfb(n, T=0, **kwargs): 

14 """GDSMFB exchange-correlation functional (spin-paired). 

15 

16 Similar to the functional with the label LDA_XC_GDSMFB and ID 577 in Libxc, but the 

17 implementations differ: https://gitlab.com/libxc/libxc/-/issues/525 

18 Exchange and correlation cannot be separated. 

19 

20 Reference: Phys. Rev. Lett. 119, 135001. 

21 

22 Args: 

23 n: Real-space electronic density. 

24 

25 Keyword Args: 

26 T: Temperature in Hartree. 

27 **kwargs: Throwaway arguments. 

28 

29 Returns: 

30 GDSMFB exchange-correlation energy density and potential. 

31 """ 

32 kwargs["zeta"] = np.zeros_like(n) 

33 exc, vxc, _ = lda_xc_gdsmfb_spin(n, T=T, **kwargs) 

34 return exc, np.array([vxc[0]]), None 

35 

36 

37def lda_xc_gdsmfb_spin(n, zeta, T=0, **kwargs): 

38 """GDSMFB exchange-correlation functional (spin-polarized). 

39 

40 Similar to the functional with the label LDA_XC_GDSMFB and ID 577 in Libxc, but the 

41 implementations differ: https://gitlab.com/libxc/libxc/-/issues/525 

42 Exchange and correlation cannot be separated. 

43 

44 Reference: Phys. Rev. Lett. 119, 135001. 

45 

46 Args: 

47 n: Real-space electronic density. 

48 zeta: Relative spin polarization. 

49 

50 Keyword Args: 

51 T: Temperature in Hartree. 

52 **kwargs: Throwaway arguments. 

53 

54 Returns: 

55 GDSMFB exchange-correlation energy density and potential. 

56 """ 

57 # Calculate properties 

58 n_up = (1 + zeta) * n / 2 

59 rs = (3 / (4 * np.pi * n)) ** (1 / 3) 

60 theta = _get_theta(T, n, zeta) 

61 theta0 = _get_theta0(theta, zeta) 

62 theta1 = _get_theta1(theta, zeta) 

63 

64 # Initialize parameters 

65 # We need to calculate coefficients for specific theta 

66 # Create a coefficients object for each theta with its corresponding parameters 

67 phi_params = PhiParams() 

68 zeta0theta = Zeta0Coeffs(theta) 

69 zeta0theta0 = Zeta0Coeffs(theta0) 

70 zeta1theta = Zeta1Coeffs(theta) 

71 zeta1theta1 = Zeta1Coeffs(theta1) 

72 

73 # Calculate fxc 

74 fxc0 = _get_fxc_zeta(rs, zeta0theta0) 

75 fxc1 = _get_fxc_zeta(rs, zeta1theta1) 

76 phi = _get_phi(rs, theta0, zeta, phi_params) 

77 fxc = fxc0 + (fxc1 - fxc0) * phi 

78 

79 # Generic derivatives 

80 drsdn = -(6 ** (1 / 3)) * (1 / n) ** (1 / 3) / (6 * np.pi ** (1 / 3) * n) 

81 dzetadn_up = -zeta / n**2 + 1 / n 

82 dzetadn_dw = -zeta / n**2 - 1 / n 

83 

84 # fxc derivatives 

85 dfxc0drs = _get_dfxc_zetadrs(rs, zeta0theta) 

86 dfxc1drs = _get_dfxc_zetadrs(rs, zeta1theta) 

87 dfxc0dtheta0 = _get_dfxc_zetadtheta(rs, zeta0theta0) 

88 dfxc1dtheta1 = _get_dfxc_zetadtheta(rs, zeta1theta1) 

89 

90 # phi derivatives 

91 dphidrs = _get_dphidrs(rs, theta0, zeta, phi_params) 

92 dphidtheta = _get_dphidtheta(rs, theta0, zeta, phi_params) 

93 dphidzeta = _get_dphidzeta(rs, theta0, zeta, phi_params) 

94 

95 # theta derivatives 

96 dthetadn_up = _get_dthetadn_up(T, n_up) 

97 dtheta0dtheta = _get_dtheta0dtheta(zeta) 

98 dtheta0dzeta = _get_dtheta0dzeta(theta0, zeta) 

99 dtheta1dtheta0 = _get_dtheta1dtheta0() 

100 

101 # Calculate vxc_up (using dndn_up=1) 

102 dfxc0a = dfxc0drs * drsdn # * dndn_up = 1 

103 dfxc1a = dfxc1drs * drsdn # * dndn_up = 1 

104 dfxc0b_up = dfxc0dtheta0 * (dtheta0dtheta * dthetadn_up + dtheta0dzeta * dzetadn_up) 

105 dfxc1b_up = ( 

106 dfxc1dtheta1 * dtheta1dtheta0 * (dtheta0dtheta * dthetadn_up + dtheta0dzeta * dzetadn_up) 

107 ) 

108 dphi_up = dphidtheta * dthetadn_up + dphidzeta * dzetadn_up + dphidrs * drsdn # * dndn_up = 1 

109 vxc_up = ( 

110 dfxc0a 

111 + dfxc0b_up 

112 - phi * (dfxc0a + dfxc0b_up - dfxc1a - dfxc1b_up) 

113 - (fxc0 - fxc1) * dphi_up 

114 ) 

115 

116 # Calculate vxc_dw (using dndn_dw=1 and dthetadn_dw=0) 

117 # dfxc0a and dfxc1a are identical for vxc_up and vxc_dw 

118 dfxc0b_dw = dfxc0dtheta0 * dtheta0dzeta * dzetadn_dw 

119 # dfxc0b += dfxc0dtheta0 * dtheta0dtheta * dthetadn_dw 

120 dfxc1b_dw = dfxc1dtheta1 * dtheta1dtheta0 * dtheta0dzeta * dzetadn_dw 

121 # dfxc1b += dfxc1dtheta1 * dtheta1dtheta0 * dtheta0dtheta * dthetadn_dw 

122 dphi_dw = dphidzeta * dzetadn_dw + dphidrs * drsdn # * dndn_dw = 1 

123 # dfxc1b += dphidtheta * dthetadn_dw 

124 vxc_dw = ( 

125 dfxc0a 

126 + dfxc0b_dw 

127 - phi * (dfxc0a + dfxc0b_dw - dfxc1a - dfxc1b_dw) 

128 - (fxc0 - fxc1) * dphi_dw 

129 ) 

130 

131 return fxc, fxc + np.array([vxc_up, vxc_dw]) * n, None 

132 

133 

134# ### Temperature dependent coefficients ### 

135 

136 

137@dataclasses.dataclass 

138class Coefficients: 

139 """Coefficients class to calculate temperature/theta dependent coefficients. 

140 

141 Reference: Phys. Rev. Lett. 119, 135001. 

142 """ 

143 

144 theta: float #: Reduced temperature. 

145 a0: float = 0.610887 

146 a1: float = 0.75 

147 a2: float = 3.04363 

148 a3: float = -0.09227 

149 a4: float = 1.7035 

150 a5: float = 8.31051 

151 a6: float = 5.1105 

152 

153 @property 

154 def b5(self): 

155 """Calculate b5.""" 

156 return self.b3 * np.sqrt(3 / 2) * self.omega * (4 / (9 * np.pi)) ** (-1 / 3) 

157 

158 @property 

159 def a(self): 

160 """Calculate a.""" 

161 with np.errstate(divide="ignore"): 

162 u = self.a0 * np.tanh( 

163 1 / self.theta, out=np.ones_like(self.theta), where=self.theta > 0 

164 ) 

165 return u * _pade(self.theta, self.a1, self.a2, self.a3, self.a4, self.a5, self.a6) 

166 

167 @property 

168 def dadtheta(self): 

169 """Calculate da / dtheta.""" 

170 with np.errstate(divide="ignore"): 

171 u = self.a0 * np.tanh( 

172 1 / self.theta, out=np.ones_like(self.theta), where=self.theta > 0 

173 ) 

174 du = np.divide( 

175 u**2 / self.a0 - self.a0, 

176 self.theta**2, 

177 out=np.zeros_like(self.theta), 

178 where=self.theta > 0, 

179 ) 

180 v = _pade(self.theta, self.a1, self.a2, self.a3, self.a4, self.a5, self.a6) 

181 dv = _dpade(self.theta, self.a1, self.a2, self.a3, self.a4, self.a5, self.a6) 

182 return du * v + u * dv 

183 

184 @property 

185 def b(self): 

186 """Calculate b.""" 

187 with np.errstate(divide="ignore"): 

188 u = np.tanh(1 / np.sqrt(self.theta), out=np.ones_like(self.theta), where=self.theta > 0) 

189 return u * _pade(self.theta, self.b1, self.b2, 0, self.b3, self.b4, self.b5) 

190 

191 @property 

192 def dbdtheta(self): 

193 """Calculate db / dtheta.""" 

194 with np.errstate(divide="ignore"): 

195 u = np.tanh(1 / np.sqrt(self.theta), out=np.ones_like(self.theta), where=self.theta > 0) 

196 du = np.divide( 

197 u**2 - 1, 

198 2 * self.theta * np.sqrt(self.theta), 

199 out=np.zeros_like(self.theta), 

200 where=self.theta > 0, 

201 ) 

202 v = _pade(self.theta, self.b1, self.b2, 0, self.b3, self.b4, self.b5) 

203 dv = _dpade(self.theta, self.b1, self.b2, 0, self.b3, self.b4, self.b5) 

204 return du * v + u * dv 

205 

206 @property 

207 def c(self): 

208 """Calculate c.""" 

209 with np.errstate(divide="ignore"): 

210 exp = np.exp(-1 / self.theta, out=np.zeros_like(self.theta), where=self.theta > 0) 

211 return (self.c1 + self.c2 * exp) * self.e 

212 

213 @property 

214 def dcdtheta(self): 

215 """Calculate dc / dtheta.""" 

216 with np.errstate(divide="ignore"): 

217 exp = np.exp(-1 / self.theta, out=np.zeros_like(self.theta), where=self.theta > 0) 

218 u = self.c1 + self.c2 * exp 

219 du = np.divide( 

220 self.c2 * exp, self.theta**2, out=np.zeros_like(self.theta), where=self.theta > 0 

221 ) 

222 v = self.e 

223 dv = self.dedtheta 

224 return du * v + u * dv 

225 

226 @property 

227 def d(self): 

228 """Calculate d.""" 

229 with np.errstate(divide="ignore"): 

230 u = np.tanh(1 / np.sqrt(self.theta), out=np.ones_like(self.theta), where=self.theta > 0) 

231 return u * _pade(self.theta, self.d1, self.d2, 0, self.d3, self.d4, self.d5) 

232 

233 @property 

234 def dddtheta(self): 

235 """Calculate dd / dtheta.""" 

236 with np.errstate(divide="ignore"): 

237 u = np.tanh(1 / np.sqrt(self.theta), out=np.ones_like(self.theta), where=self.theta > 0) 

238 du = np.divide( 

239 u**2 - 1, 

240 2 * self.theta * np.sqrt(self.theta), 

241 out=np.zeros_like(self.theta), 

242 where=self.theta > 0, 

243 ) 

244 v = _pade(self.theta, self.d1, self.d2, 0, self.d3, self.d4, self.d5) 

245 dv = _dpade(self.theta, self.d1, self.d2, 0, self.d3, self.d4, self.d5) 

246 return du * v + u * dv 

247 

248 @property 

249 def e(self): 

250 """Calculate e.""" 

251 with np.errstate(divide="ignore"): 

252 u = np.tanh(1 / self.theta, out=np.ones_like(self.theta), where=self.theta > 0) 

253 return u * _pade(self.theta, self.e1, self.e2, 0, self.e3, self.e4, self.e5) 

254 

255 @property 

256 def dedtheta(self): 

257 """Calculate de / dtheta.""" 

258 with np.errstate(divide="ignore"): 

259 u = np.tanh(1 / self.theta, out=np.ones_like(self.theta), where=self.theta > 0) 

260 du = np.divide(u**2 - 1, self.theta**2, out=np.zeros_like(self.theta), where=self.theta > 0) 

261 v = _pade(self.theta, self.e1, self.e2, 0, self.e3, self.e4, self.e5) 

262 dv = _dpade(self.theta, self.e1, self.e2, 0, self.e3, self.e4, self.e5) 

263 return du * v + u * dv 

264 

265 

266# ### Parameters ### 

267 

268 

269@dataclasses.dataclass 

270class Zeta0Coeffs(Coefficients): 

271 """Coefficient class using the parameters for zeta=0. 

272 

273 Reference: Phys. Rev. Lett. 119, 135001. 

274 """ 

275 

276 omega: float = 1 

277 b1: float = 0.3436902 

278 b2: float = 7.82159531356 

279 b3: float = 0.300483986662 

280 b4: float = 15.8443467125 

281 c1: float = 0.8759442 

282 c2: float = -0.230130843551 

283 d1: float = 0.72700876 

284 d2: float = 2.38264734144 

285 d3: float = 0.30221237251 

286 d4: float = 4.39347718395 

287 d5: float = 0.729951339845 

288 e1: float = 0.25388214 

289 e2: float = 0.815795138599 

290 e3: float = 0.0646844410481 

291 e4: float = 15.0984620477 

292 e5: float = 0.230761357474 

293 

294 

295@dataclasses.dataclass 

296class Zeta1Coeffs(Coefficients): 

297 """Coefficient class using the parameters for zeta=1. 

298 

299 Reference: Phys. Rev. Lett. 119, 135001. 

300 """ 

301 

302 omega: float = 2 ** (1 / 3) 

303 b1: float = 0.84987704 

304 b2: float = 3.04033012073 

305 b3: float = 0.0775730131248 

306 b4: float = 7.57703592489 

307 c1: float = 0.91126873 

308 c2: float = -0.0307957123308 

309 d1: float = 1.48658718 

310 d2: float = 4.92684905511 

311 d3: float = 0.0849387225179 

312 d4: float = 8.3269821188 

313 d5: float = 0.218864952126 

314 e1: float = 0.27454097 

315 e2: float = 0.400994856555 

316 e3: float = 2.88773194962 

317 e4: float = 6.33499237092 

318 e5: float = 24.823008753 

319 

320 

321@dataclasses.dataclass 

322class PhiParams: 

323 """Parameter class holding the spin-interpolation function parameters. 

324 

325 Reference: Phys. Rev. Lett. 119, 135001. 

326 """ 

327 

328 # Sign of parameters is different from the supplemental material 

329 h1: float = 3.18747258 

330 h2: float = 7.74662802 

331 lambda1: float = 1.85909536 

332 lambda2: float = 0 

333 

334 

335# ### Pade approximation and derivative ### 

336 

337 

338def _pade(x, n1, n2, n3, n4, d1, d2): 

339 """Pade approximation. 

340 

341 Not the general case but as needed for this functional. 

342 """ 

343 num = n1 + n2 * x**2 + n3 * x**3 + n4 * x**4 

344 denom = 1 + d1 * x**2 + d2 * x**4 

345 return num / denom 

346 

347 

348def _dpade(x, n1, n2, n3, n4, d1, d2): 

349 """Pade approximation derivative.""" 

350 num = n1 + n2 * x**2 + n3 * x**3 + n4 * x**4 

351 denom = 1 + d1 * x**2 + d2 * x**4 

352 dnum = 2 * n2 * x + 3 * n3 * x**2 + 4 * n4 * x**3 

353 ddenom = 2 * d1 * x + 4 * d2 * x**3 

354 # df = (a'b - ab') / b^2 

355 return (dnum * denom - num * ddenom) / denom**2 

356 

357 

358# ### theta and derivatives ### 

359 

360 

361def _get_theta(T, n, zeta): 

362 """Calculate the reduced temperature theta. 

363 

364 Reference: Phys. Rev. Lett. 119, 135001. 

365 Only mentioned in the arXiv version: https://arxiv.org/abs/1703.08074 

366 """ 

367 n_up = (1 + zeta) * n / 2 

368 T_fermi = (6 * np.pi**2 * n_up) ** (2 / 3) / 2 

369 return T / T_fermi 

370 

371 

372def _get_dthetadn_up(T, n_up): 

373 """Calculate dtheta / dn_up.""" 

374 return -4 / (3 * 6 ** (2 / 3)) * T / (np.pi ** (4 / 3) * n_up ** (5 / 3)) 

375 

376 

377def _get_theta0(theta, zeta): 

378 """Calculate theta0. 

379 

380 Reference: Phys. Rev. Lett. 119, 135001. 

381 """ 

382 return theta * (1 + zeta) ** (2 / 3) 

383 

384 

385def _get_dtheta0dtheta(zeta): 

386 """Calculate dtheta0 / dtheta.""" 

387 return (1 + zeta) ** (2 / 3) 

388 

389 

390def _get_dtheta0dzeta(theta, zeta): 

391 """Calculate dtheta0 / dzeta.""" 

392 return 2 * theta / (3 * (1 + zeta) ** (1 / 3)) 

393 

394 

395def _get_theta1(theta, zeta): 

396 """Calculate theta1. 

397 

398 Reference: Phys. Rev. Lett. 119, 135001. 

399 

400 It is not explicitly mentioned but introduced as used in Eq. (5). 

401 """ 

402 theta0 = _get_theta0(theta, zeta) 

403 return theta0 * 2 ** (-2 / 3) 

404 

405 

406def _get_dtheta1dtheta0(): 

407 """Calculate dtheta1 / dtheta0.""" 

408 return 2 ** (-2 / 3) 

409 

410 

411# ### fxc_zeta and derivatives ### 

412 

413 

414def _get_fxc_zeta(rs, p): 

415 """Calculate the Pade formula f_xc^zeta. 

416 

417 Reference: Phys. Rev. Lett. 119, 135001. 

418 """ 

419 num = p.omega * p.a + np.sqrt(rs) * p.b + rs * p.c 

420 denom = 1 + np.sqrt(rs) * p.d + rs * p.e 

421 return -1 / rs * num / denom 

422 

423 

424def _get_dfxc_zetadrs(rs, p): 

425 """Calculate dfxc_zeta / drs.""" 

426 tmp1 = (-p.b / (2 * np.sqrt(rs)) - p.c) / (rs * (p.d * np.sqrt(rs) + p.e * rs + 1)) 

427 tmp2 = (-p.d / (2 * np.sqrt(rs)) - p.e) * (-p.a * p.omega - p.b * np.sqrt(rs) - p.c * rs) 

428 tmp3 = rs * (p.d * np.sqrt(rs) + p.e * rs + 1) ** 2 

429 tmp4 = (-p.a * p.omega - p.b * np.sqrt(rs) - p.c * rs) / ( 

430 rs**2 * (p.d * np.sqrt(rs) + p.e * rs + 1) 

431 ) 

432 return tmp1 + tmp2 / tmp3 - tmp4 

433 

434 

435def _get_dfxc_zetadtheta(rs, p): 

436 """Calculate dfxc_zeta / dzeta.""" 

437 tmp1 = (-np.sqrt(rs) * p.dddtheta - rs * p.dedtheta) * ( 

438 -p.omega * p.a - p.b * np.sqrt(rs) - p.c * rs 

439 ) 

440 tmp2 = (p.d * np.sqrt(rs) + p.e * rs + 1) ** 2 * rs 

441 tmp3 = -p.omega * p.dadtheta - np.sqrt(rs) * p.dbdtheta - rs * p.dcdtheta 

442 tmp4 = (p.d * np.sqrt(rs) + p.e * rs + 1) * rs 

443 return tmp1 / tmp2 + tmp3 / tmp4 

444 

445 

446# ### phi and derivatives ### 

447 

448 

449def _get_phi(rs, theta, zeta, phi_params): 

450 """Calculate the interpolation function phi. 

451 

452 Reference: Phys. Rev. Lett. 119, 135001. 

453 """ 

454 alpha = _get_alpha(rs, theta, phi_params) 

455 return ((1 + zeta) ** alpha + (1 - zeta) ** alpha - 2) / (2**alpha - 2) 

456 

457 

458def _get_dphidrs(rs, theta, zeta, phi_params): 

459 """Calculate dphi / drs.""" 

460 alpha = _get_alpha(rs, theta, phi_params) 

461 tmp1 = (1 - zeta) ** alpha * np.log(1 - zeta, out=np.zeros_like(zeta), where=1 - zeta > 0) 

462 tmp2 = (1 + zeta) ** alpha * np.log(1 + zeta) 

463 duv = (tmp1 + tmp2) * (2**alpha - 2) 

464 udv = ((1 - zeta) ** alpha + (1 + zeta) ** alpha - 2) * 2**alpha * np.log(2) 

465 vv = (2**alpha - 2) ** 2 

466 dalphadrs = _get_dalphadrs(rs, theta, phi_params) 

467 return (duv - udv) * dalphadrs / vv 

468 

469 

470def _get_dphidtheta(rs, theta, zeta, phi_params): 

471 """Calculate dphi / dtheta.""" 

472 alpha = _get_alpha(rs, theta, phi_params) 

473 dalphadtheta = _get_dalphadtheta(rs, theta, phi_params) 

474 tmp1 = (1 - zeta) ** alpha * np.log(1 - zeta, out=np.zeros_like(zeta), where=1 - zeta > 0) 

475 tmp2 = (1 + zeta) ** alpha * np.log(1 + zeta) 

476 duv = (tmp1 + tmp2) * (2**alpha - 2) 

477 udv = ((1 - zeta) ** alpha + (1 + zeta) ** alpha - 2) * 2**alpha * np.log(2) 

478 vv = (2**alpha - 2) ** 2 

479 dalphadtheta = _get_dalphadtheta(rs, theta, phi_params) 

480 return (duv - udv) * dalphadtheta / vv 

481 

482 

483def _get_dphidzeta(rs, theta, zeta, phi_params): 

484 """Calculate dphi / dzeta.""" 

485 alpha = _get_alpha(rs, theta, phi_params) 

486 tmp1 = alpha * (1 + zeta) ** alpha / (1 + zeta) 

487 tmp2 = np.divide( 

488 alpha * (1 - zeta) ** alpha, 1 - zeta, out=np.zeros_like(zeta), where=1 - zeta > 0 

489 ) 

490 return (tmp1 - tmp2) / (2**alpha - 2) 

491 

492 

493# ### alpha and derivatives ### 

494 

495 

496def _get_alpha(rs, theta, phi_params): 

497 """Calculate alpha. 

498 

499 Reference: Phys. Rev. Lett. 119, 135001. 

500 """ 

501 h = _get_h(rs, phi_params) 

502 lamda = _get_lambda(rs, theta, phi_params) 

503 return 2 - h * np.exp(-theta * lamda) 

504 

505 

506def _get_dalphadrs(rs, theta, phi_params): 

507 """Calculate dalpha / drs.""" 

508 h = _get_h(rs, phi_params) 

509 lamda = _get_lambda(rs, theta, phi_params) 

510 dhdrs = _get_dhdrs(rs, phi_params) 

511 dlambdadrs = _get_dlambdadrs(rs, theta, phi_params) 

512 return -dhdrs * np.exp(-theta * lamda) + dlambdadrs * theta * h * np.exp(-theta * lamda) 

513 

514 

515def _get_dalphadtheta(rs, theta, phi_params): 

516 """Calculate dalpha / dtheta.""" 

517 h = _get_h(rs, phi_params) 

518 lamda = _get_lambda(rs, theta, phi_params) 

519 dlambdadtheta = _get_dlambdadtheta(rs, phi_params) 

520 return (dlambdadtheta * theta + lamda) * h * np.exp(-theta * lamda) 

521 

522 

523# ### h and derivative ### 

524 

525 

526def _get_h(rs, phi_params): 

527 """Calculate h. 

528 

529 Reference: Phys. Rev. Lett. 119, 135001. 

530 """ 

531 return (2 / 3 + phi_params.h1 * rs) / (1 + phi_params.h2 * rs) 

532 

533 

534def _get_dhdrs(rs, phi_params): 

535 """Calculate dh / drs.""" 

536 return ( 

537 phi_params.h1 / (phi_params.h2 * rs + 1) 

538 - phi_params.h2 * (phi_params.h1 * rs + 2 / 3) / (phi_params.h2 * rs + 1) ** 2 

539 ) 

540 

541 

542# ### lambda and derivatives ### 

543 

544 

545def _get_lambda(rs, theta, phi_params): 

546 """Calculate lambda. 

547 

548 Reference: Phys. Rev. Lett. 119, 135001. 

549 """ 

550 return phi_params.lambda1 + phi_params.lambda2 * theta * rs ** (1 / 2) 

551 

552 

553def _get_dlambdadrs(rs, theta, phi_params): 

554 """Calculate dlambda / drs.""" 

555 return phi_params.lambda2 * theta / (2 * np.sqrt(rs)) 

556 

557 

558def _get_dlambdadtheta(rs, phi_params): 

559 """Calculate dlambda / dtheta.""" 

560 return phi_params.lambda2 * np.sqrt(rs)