Coverage for eminus/tools.py: 96.12%

1# SPDX-FileCopyrightText: 2021 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""Various tools to check physical properties.""" 

4 

5import math 

6import numbers 

7 

8import numpy as np 

9from scipy.optimize import minimize_scalar, root_scalar 

10 

11from . import backend as xp 

12from .dft import get_epsilon, get_epsilon_unocc 

13from .gga import get_grad_field, get_tau 

14from .logger import log 

15from .utils import handle_k 

16 

17 

18def cutoff2gridspacing(E): 

19 """Convert plane wave energy cut-off to a real-space grid spacing. 

20 

21 Reference: Phys. Rev. B 54, 14362. 

22 

23 Args: 

24 E: Energy in Hartree. 

25 

26 Returns: 

27 Grid spacing in Bohr. 

28 """ 

29 return math.pi / math.sqrt(2 * E) 

30 

31 

32def gridspacing2cutoff(h): 

33 """Convert real-space grid spacing to plane wave energy cut-off. 

34 

35 Reference: Phys. Rev. B 54, 14362. 

36 

37 Args: 

38 h: Grid spacing in Bohr. 

39 

40 Returns: 

41 Cut-off in Hartree. 

42 """ 

43 return 0.5 * (math.pi / h) ** 2 

44 

45 

46def center_of_mass(coords, masses=None): 

47 """Calculate the center of mass for a set of coordinates and masses. 

48 

49 Args: 

50 coords: Array of real-space coordinates. 

51 

52 Keyword Args: 

53 masses: Mass or weight for each coordinate. 

54 

55 Returns: 

56 Center of mass. 

57 """ 

58 if coords is None: 

59 msg = 'The provided coordinates are "None".' 

60 raise ValueError(msg) 

61 if masses is None: 

62 masses = xp.ones(len(coords)) 

63 return xp.sum(masses * coords.T, axis=1) / xp.sum(masses) 

64 

65 

66@handle_k 

67def orbital_center(obj, psirs): 

68 """Calculate the orbital center of masses, e.g., from localized orbitals. 

69 

70 Args: 

71 obj: Atoms or SCF object. 

72 psirs: Set of orbitals in real-space. 

73 

74 Returns: 

75 Center of masses. 

76 """ 

77 atoms = obj._atoms 

78 

79 coms = [xp.empty((0, 3))] * 2 

80 Ncom = psirs.shape[2] 

81 for spin in range(atoms.occ.Nspin): 

82 coms_spin = xp.empty((Ncom, 3)) 

83 

84 # Squared orbitals 

85 psi2 = xp.real(psirs[spin].conj() * psirs[spin]) 

86 for i in range(Ncom): 

87 coms_spin[i] = center_of_mass(atoms.r, psi2[:, i]) 

88 coms[spin] = coms_spin 

89 return coms 

90 

91 

92def inertia_tensor(coords, masses=None): 

93 """Calculate the inertia tensor for a set of coordinates and masses. 

94 

95 Reference: https://en.wikipedia.org/wiki/Moment_of_inertia 

96 

97 Args: 

98 coords: Array of real-space coordinates. 

99 

100 Keyword Args: 

101 masses: Mass or weight for each coordinate. 

102 

103 Returns: 

104 Inertia tensor. 

105 """ 

106 if masses is None: 

107 masses = xp.ones(len(coords)) 

108 

109 # The inertia tensor for a set of point masses can be calculated with a simple summation 

110 # https://en.wikipedia.org/wiki/Moment_of_inertia#Definition_2 

111 I = xp.empty((3, 3)) 

112 I[0, 0] = xp.sum(masses * (coords[:, 1] ** 2 + coords[:, 2] ** 2)) 

113 I[1, 1] = xp.sum(masses * (coords[:, 0] ** 2 + coords[:, 2] ** 2)) 

114 I[2, 2] = xp.sum(masses * (coords[:, 0] ** 2 + coords[:, 1] ** 2)) 

115 I[0, 1] = I[1, 0] = -xp.sum(masses * (coords[:, 0] * coords[:, 1])) 

116 I[0, 2] = I[2, 0] = -xp.sum(masses * (coords[:, 0] * coords[:, 2])) 

117 I[1, 2] = I[2, 1] = -xp.sum(masses * (coords[:, 1] * coords[:, 2])) 

118 return I 

119 

120 

121def get_dipole(scf, n=None): 

122 """Calculate the electric dipole moment. 

123 

124 This function does not account for periodicity, it may be a good idea to center the system. 

125 

126 Reference: J. Chem. Phys. 155, 224109. 

127 

128 Args: 

129 scf: SCF object. 

130 

131 Keyword Args: 

132 n: Real-space electronic density. 

133 

134 Returns: 

135 Electric dipole moment in e * Bohr. 

136 """ 

137 atoms = scf.atoms 

138 if n is None: 

139 if scf.n is None: 

140 log.error("There is no density to calculate a dipole moment.") 

141 return 0 

142 n = scf.n 

143 

144 # Diple moment: mu = \sum Z pos - \int n(r) r dr 

145 mu = xp.zeros(3) 

146 for i in range(atoms.Natoms): 

147 mu += atoms.Z[i] * atoms.pos[i] 

148 

149 for dim in range(3): 

150 mu[dim] -= atoms.dV * xp.sum(n * atoms.r[:, dim]) 

151 return mu 

152 

153 

154def get_ip(scf): 

155 """Calculate the ionization potential by calculating the negative HOMO energy. 

156 

157 Reference: Physica 1, 104. 

158 

159 Args: 

160 scf: SCF object. 

161 

162 Returns: 

163 Ionization potential in Hartree. 

164 """ 

165 scf.kpts._assert_gamma_only() 

166 epsilon = get_epsilon(scf, scf.W)[0] 

167 # Account for spin-polarized calculations 

168 epsilon = xp.sort(xp.ravel(epsilon)) 

169 return -epsilon[-1] 

170 

171 

172@handle_k(mode="skip") 

173def check_ortho(obj, func, eps=1e-9): 

174 """Check the orthogonality condition for a set of functions. 

175 

176 Args: 

177 obj: Atoms or SCF object. 

178 func: A discretized set of functions. 

179 

180 Keyword Args: 

181 eps: Tolerance for the condition. 

182 

183 Returns: 

184 Orthogonality status for the set of functions. 

185 """ 

186 atoms = obj._atoms 

187 func = xp.atleast_3d(func) 

188 

189 # It makes no sense to calculate anything for only one function 

190 if atoms.occ.Nstate == 1: 

191 log.warning("Need at least two functions to check their orthogonality.") 

192 return True 

193 

194 ortho_bool = True 

195 # Check the condition for every combination 

196 # Orthogonality condition: \int func1^* func2 dr = 0 

197 for spin in range(atoms.occ.Nspin): 

198 for i in range(atoms.occ.Nstate): 

199 for j in range(i + 1, atoms.occ.Nstate): 

200 res = atoms.dV * xp.sum(func[spin, :, i].conj() * func[spin, :, j]) 

201 tmp_bool = abs(res) < eps 

202 ortho_bool *= tmp_bool 

203 log.debug(f"Function {i} and {j}:\nValue: {res:.7f}\nOrthogonal: {tmp_bool}") 

204 log.info(f"Orthogonal: {ortho_bool}") 

205 return ortho_bool 

206 

207 

208@handle_k(mode="skip") 

209def check_norm(obj, func, eps=1e-9): 

210 """Check the normalization condition for a set of functions. 

211 

212 Args: 

213 obj: Atoms or SCF object. 

214 func: A discretized set of functions. 

215 

216 Keyword Args: 

217 eps: Tolerance for the condition. 

218 

219 Returns: 

220 Normalization status for the set of functions. 

221 """ 

222 atoms = obj._atoms 

223 func = xp.atleast_3d(func) 

224 

225 norm_bool = True 

226 # Check the condition for every function 

227 # Normality condition: \int func^* func dr = 1 

228 for spin in range(atoms.occ.Nspin): 

229 for i in range(atoms.occ.Nstate): 

230 res = atoms.dV * xp.sum(func[spin, :, i].conj() * func[spin, :, i]) 

231 tmp_bool = abs(1 - res) < eps 

232 norm_bool *= tmp_bool 

233 log.debug(f"Function {i}:\nValue: {res:.7f}\nNormalized: {tmp_bool}") 

234 log.info(f"Normalized: {norm_bool}") 

235 return norm_bool 

236 

237 

238@handle_k(mode="skip") 

239def check_orthonorm(obj, func, eps=1e-9): 

240 """Check the orthonormality conditions for a set of functions. 

241 

242 Args: 

243 obj: Atoms or SCF object. 

244 func: A discretized set of functions. 

245 

246 Keyword Args: 

247 eps: Tolerance for the condition. 

248 

249 Returns: 

250 Orthonormality status for the set of functions. 

251 """ 

252 atoms = obj._atoms 

253 ortho_bool = check_ortho(atoms, func, eps) 

254 norm_bool = check_norm(atoms, func, eps) 

255 log.info(f"Orthonormal: {ortho_bool * norm_bool}") 

256 return ortho_bool * norm_bool 

257 

258 

259def get_isovalue(n, percent=85): 

260 """Find an isovalue that contains a percentage of the electronic density. 

261 

262 Reference: J. Chem. Phys. 158, 164102. 

263 

264 Args: 

265 n: Real-space electronic density. 

266 

267 Keyword Args: 

268 percent: Amount of density that should be contained. 

269 

270 Returns: 

271 Isovalue that contains the specified percentage of the density. 

272 """ 

273 

274 def deviation(isovalue): 

275 """Wrapper function for finding the isovalue by minimization.""" 

276 n_mask = xp.sum(n[n > isovalue]) 

277 return float(xp.abs(percent - (n_mask / n_ref) * 100)) 

278 

279 # Integrated density 

280 n_ref = xp.sum(n) 

281 # Finding the isovalue is an optimization problem, minimizing the deviation above 

282 # The problem is bound by zero (no density) and the maximum value in n 

283 res = minimize_scalar(deviation, bounds=(0, float(xp.max(n))), method="bounded") 

284 return res.x 

285 

286 

287def get_tautf(scf): 

288 """Calculate the Thomas-Fermi kinetic energy densities per spin. 

289 

290 Reference: Phys. Lett. B 63, 395. 

291 

292 Args: 

293 scf: SCF object. 

294 

295 Returns: 

296 Real-space Thomas-Fermi kinetic energy density. 

297 """ 

298 atoms = scf.atoms 

299 # Use the definition with a division by two 

300 tautf = 3 / 10 * (atoms.occ.Nspin * 3 * math.pi**2) ** (2 / 3) * scf.n_spin ** (5 / 3) 

301 

302 log.debug(f"Calculated Ekin: {scf.energies.Ekin:.6f} Eh") 

303 log.debug(f"Integrated tautf: {xp.sum(tautf) * atoms.dV:.6f} Eh") 

304 return tautf 

305 

306 

307def get_tauw(scf): 

308 """Calculate the von Weizsaecker kinetic energy densities per spin. 

309 

310 Reference: Z. Phys. 96, 431. 

311 

312 Args: 

313 scf: SCF object. 

314 

315 Returns: 

316 Real-space von Weizsaecker kinetic energy density. 

317 """ 

318 atoms = scf.atoms 

319 if scf.dn_spin is None: 

320 dn_spin = get_grad_field(atoms, scf.n_spin) 

321 else: 

322 dn_spin = scf.dn_spin 

323 dn2 = xp.linalg.norm(dn_spin, axis=2) ** 2 

324 # Use the definition with a division by two 

325 tauw = dn2 / (8 * scf.n_spin) 

326 

327 # For one- and two-electron systems the integrated KED has to be the same as the calculated KE 

328 log.debug(f"Calculated Ekin: {scf.energies.Ekin:.6f} Eh") 

329 log.debug(f"Integrated tauw: {xp.sum(tauw) * atoms.dV:.6f} Eh") 

330 return tauw 

331 

332 

333def get_elf(scf): 

334 """Calculate the electron localization function. 

335 

336 Reference: J. Chem. Phys. 92, 5397. 

337 

338 Args: 

339 scf: SCF object. 

340 

341 Returns: 

342 Real-space electron localization function. 

343 """ 

344 D = get_tau(scf.atoms, scf.Y) - get_tauw(scf) 

345 D0 = get_tautf(scf) 

346 X = D / D0 

347 return 1 / (1 + X**2) 

348 

349 

350def get_reduced_gradient(scf, eps=0): 

351 """Calculate the reduced density gradient s. 

352 

353 Reference: Phys. Rev. Lett. 77, 3865. 

354 

355 Args: 

356 scf: SCF object. 

357 

358 Kwargs: 

359 eps: Threshold of the density where s should be truncated. 

360 

361 Returns: 

362 Real-space reduced density gradient. 

363 """ 

364 atoms = scf.atoms 

365 if scf.dn_spin is None: 

366 dn_spin = get_grad_field(atoms, scf.n_spin) 

367 else: 

368 dn_spin = scf.dn_spin 

369 norm_dn = xp.linalg.norm(xp.sum(dn_spin, axis=0), axis=1) 

370 

371 kf = (3 * math.pi**2 * scf.n) ** (1 / 3) 

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

373 s = norm_dn / (2 * kf * scf.n) 

374 s[scf.n < eps] = 0 

375 return s 

376 

377 

378def get_spin_squared(scf): 

379 """Calculate the expectation value of the squared spin operator <S^2>. 

380 

381 Reference: Appl. Phys. Express 12, 115506. 

382 

383 Args: 

384 scf: SCF object. 

385 

386 Returns: 

387 The DFT value for <S^2>. 

388 """ 

389 atoms = scf.atoms 

390 # <S^2> for a restricted calculation is always zero 

391 if not atoms.unrestricted: 

392 return 0 

393 

394 rhoXr = scf.n_spin[0] - scf.n_spin[1] 

395 rhoXr[rhoXr < 0] = 0 

396 rhoX = xp.sum(rhoXr) * atoms.dV 

397 SX = 0.5 * (xp.sum(scf.n_spin[0]) - xp.sum(scf.n_spin[1])) * atoms.dV 

398 return SX * (SX + 1) + rhoX 

399 

400 

401def get_multiplicity(scf): 

402 """Calculate the multiplicity from <S^2>. 

403 

404 Args: 

405 scf: SCF object. 

406 

407 Returns: 

408 Multiplicity 2S+1. 

409 """ 

410 S2 = get_spin_squared(scf) 

411 # <S^2> = S(S+1) = S^2+S+0.25-0.25 = (S+0.5)^2-0.25 => S = sqrt(<S^2>+0.25)-0.5 

412 S = math.sqrt(S2 + 0.25) - 0.5 

413 return 2 * S + 1 

414 

415 

416def get_magnetization(scf): 

417 """Calculate the total magnetization M. 

418 

419 Args: 

420 scf: SCF object. 

421 

422 Returns: 

423 Total magnetization. 

424 """ 

425 # For a spin paired calculation the total magnetization is just zero 

426 if not scf.atoms.unrestricted: 

427 return 0 

428 

429 return xp.sum(scf.n_spin[0] - scf.n_spin[1]) / xp.sum(scf.n) 

430 

431 

432def get_bandgap(scf): 

433 """Calculate the band gap. 

434 

435 Args: 

436 scf: SCF object. 

437 

438 Returns: 

439 Band gap energy. 

440 """ 

441 e_occ = get_epsilon(scf, scf.W, **scf._precomputed) 

442 

443 if scf.Z is None: 

444 log.warning("The SCF object has no unoccupied energies, can't calculate band gap.") 

445 return 0 

446 

447 e_unocc = get_epsilon_unocc(scf, scf.W, scf.Z, **scf._precomputed) 

448 return xp.min(e_unocc) - xp.max(e_occ) 

449 

450 

451def get_Efermi(obj, epsilon=None): 

452 """Calculate the Fermi energy. 

453 

454 Reference: Phys. Rev. B 107, 195122. 

455 

456 Args: 

457 obj: SCF or Occupations object. 

458 

459 Keyword Args: 

460 epsilon: Eigenenergies. 

461 

462 Returns: 

463 Fermi energy. 

464 """ 

465 # Handle the obj argument 

466 if hasattr(obj, "smearing"): 

467 if epsilon is None: 

468 log.error("When passing an Occupations object the eigenenergies have to be given.") 

469 occ = obj 

470 else: 

471 occ = obj.atoms.occ 

472 

473 # Calculate the eigenenergies if necessary 

474 if epsilon is None: 

475 e_occ = get_epsilon(obj, obj.W, **obj._precomputed) 

476 else: 

477 e_occ = epsilon 

478 

479 def electron_root(Efermi): 

480 """Number of electrons by Fermi distribution minus the actual number of electrons.""" 

481 occ_sum = 0 

482 for ik in range(occ.Nk): 

483 occ_sum += occ.wk[ik] * xp.sum(fermi_distribution(e_occ[ik], Efermi, occ.smearing)) 

484 return float(occ_sum * 2 / occ.Nspin - occ.Nelec) 

485 

486 # For smeared systems we have to find the root of an objective function 

487 if occ.smearing > 0: 

488 return root_scalar(electron_root, bracket=(float(xp.min(e_occ)), float(xp.max(e_occ)))).root 

489 

490 if obj.Z is None: 

491 log.warning("The SCF object has no unoccupied energies, return the maximum energy instead.") 

492 return xp.max(e_occ) 

493 

494 e_unocc = get_epsilon_unocc(obj, obj.W, obj.Z, **obj._precomputed) 

495 return xp.max(e_occ) + (xp.min(e_unocc) - xp.max(e_occ)) / 2 

496 

497 

498def fermi_distribution(E, mu, kbT): 

499 """Calculate the Fermi distribution. 

500 

501 Reference: https://en.wikipedia.org/wiki/Fermi%E2%80%93Dirac_statistics 

502 

503 Args: 

504 E: State energy. 

505 mu: Chemical energy or Fermi energy. 

506 kbT: Thermic energy or smearing width. 

507 

508 Returns: 

509 Fermi distribution. 

510 """ 

511 x = (E - mu) / kbT 

512 with np.errstate(over="ignore"): 

513 if isinstance(x, numbers.Real): 

514 return 1 / (math.exp(x) + 1) 

515 return 1 / (xp.exp(x) + 1) 

516 

517 

518def electronic_entropy(E, mu, kbT): 

519 """Calculate the electronic entropic energy. 

520 

521 Reference: J. Phys. Condens. Matter 1, 689. 

522 

523 Args: 

524 E: State energy. 

525 mu: Chemical energy or Fermi energy. 

526 kbT: Thermic energy or smearing width. 

527 

528 Returns: 

529 Electronic entropic energy. 

530 """ 

531 # Condition taken from: https://gitlab.com/QEF/q-e/-/blob/master/Modules/w1gauss.f90 

532 if abs((E - mu) / kbT) > 36: 

533 return 0 

534 f = fermi_distribution(E, mu, kbT) 

535 if isinstance(f, numbers.Real): 

536 return f * math.log(f) + (1 - f) * math.log(1 - f) 

537 return f * xp.log(f) + (1 - f) * xp.log(1 - f) 

538 

539 

540def get_dos(epsilon, wk, spin=0, npts=500, width=0.1): 

541 """Calculate the total density of states. 

542 

543 Reference: https://gitlab.com/gpaw/gpaw/-/blob/master/gpaw/calculator.py 

544 

545 Args: 

546 epsilon: Eigenenergies. 

547 wk: Chemical energy or Fermi energy. 

548 

549 Keyword Args: 

550 spin: Spin channel. 

551 npts: Number of energy discretizations. 

552 width: Gaussian width. 

553 

554 Returns: 

555 Eigenenergies and DOS. 

556 """ 

557 

558 def delta(x, x0, width): 

559 """Gaussian of given width centered at x0.""" 

560 return xp.exp(-(((x - x0) / width) ** 2)) / (math.sqrt(math.pi) * width) 

561 

562 energies = epsilon[:, spin].flatten() 

563 emin = xp.min(energies) - 5 * width 

564 emax = xp.max(energies) + 5 * width 

565 e = xp.linspace(emin, emax, npts) 

566 dos_e = xp.zeros(npts) 

567 for e0, w in zip(energies, wk): 

568 dos_e += w * delta(e, e0, width) 

569 return e, dos_e