Coverage for eminus/tools.py: 96.50%

1# SPDX-FileCopyrightText: 2021 The eminus developers

2# SPDX-License-Identifier: Apache-2.0

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

5import numpy as np

6from scipy.linalg import norm

7from scipy.optimize import minimize_scalar, root_scalar

9from .dft import get_epsilon, get_epsilon_unocc

10from .gga import get_grad_field, get_tau

11from .logger import log

12from .utils import handle_k

15def cutoff2gridspacing(E):

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

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

20 Args:

21 E: Energy in Hartree.

23 Returns:

24 Grid spacing in Bohr.

25 """

26 return np.pi / np.sqrt(2 * E)

29def gridspacing2cutoff(h):

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

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

34 Args:

35 h: Grid spacing in Bohr.

37 Returns:

38 Cut-off in Hartree.

39 """

40 return 0.5 * (np.pi / h) ** 2

43def center_of_mass(coords, masses=None):

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

46 Args:

47 coords: Array of real-space coordinates.

49 Keyword Args:

50 masses: Mass or weight for each coordinate.

52 Returns:

53 Center of mass.

54 """

55 if coords is None:

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

57 raise ValueError(msg)

58 if masses is None:

59 masses = np.ones(len(coords))

60 return np.sum(masses * coords.T, axis=1) / np.sum(masses)

63@handle_k

64def orbital_center(obj, psirs):

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

67 Args:

68 obj: Atoms or SCF object.

69 psirs: Set of orbitals in real-space.

71 Returns:

72 Center of masses.

73 """

74 atoms = obj._atoms

76 coms = [np.empty((0, 3))] * 2

77 Ncom = psirs.shape[2]

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

79 coms_spin = np.empty((Ncom, 3))

81 # Squared orbitals

82 psi2 = np.real(psirs[spin].conj() * psirs[spin])

83 for i in range(Ncom):

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

85 coms[spin] = coms_spin

86 return coms

89def inertia_tensor(coords, masses=None):

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

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

94 Args:

95 coords: Array of real-space coordinates.

97 Keyword Args:

98 masses: Mass or weight for each coordinate.

100 Returns:

101 Inertia tensor.

102 """

103 if masses is None:

104 masses = np.ones(len(coords))

105

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

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

108 I = np.empty((3, 3))

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

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

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

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

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

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

115 return I

116

117

118def get_dipole(scf, n=None):

119 """Calculate the electric dipole moment.

120

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

122

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

124

125 Args:

126 scf: SCF object.

127

128 Keyword Args:

129 n: Real-space electronic density.

130

131 Returns:

132 Electric dipole moment in e * Bohr.

133 """

134 atoms = scf.atoms

135 if n is None:

136 if scf.n is None:

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

138 return 0

139 n = scf.n

140

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

142 mu = np.zeros(3)

143 for i in range(atoms.Natoms):

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

145

146 for dim in range(3):

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

148 return mu

149

150

151def get_ip(scf):

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

153

154 Reference: Physica 1, 104.

155

156 Args:

157 scf: SCF object.

158

159 Returns:

160 Ionization potential in Hartree.

161 """

162 scf.kpts._assert_gamma_only()

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

164 # Account for spin-polarized calculations

165 epsilon = np.sort(np.ravel(epsilon))

166 return -epsilon[-1]

167

168

169@handle_k(mode="skip")

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

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

172

173 Args:

174 obj: Atoms or SCF object.

175 func: A discretized set of functions.

176

177 Keyword Args:

178 eps: Tolerance for the condition.

179

180 Returns:

181 Orthogonality status for the set of functions.

182 """

183 atoms = obj._atoms

184 func = np.atleast_3d(func)

185

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

187 if atoms.occ.Nstate == 1:

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

189 return True

190

191 ortho_bool = True

192 # Check the condition for every combination

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

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

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

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

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

198 tmp_bool = abs(res) < eps

199 ortho_bool *= tmp_bool

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

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

202 return ortho_bool

203

204

205@handle_k(mode="skip")

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

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

208

209 Args:

210 obj: Atoms or SCF object.

211 func: A discretized set of functions.

212

213 Keyword Args:

214 eps: Tolerance for the condition.

215

216 Returns:

217 Normalization status for the set of functions.

218 """

219 atoms = obj._atoms

220 func = np.atleast_3d(func)

221

222 norm_bool = True

223 # Check the condition for every function

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

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

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

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

228 tmp_bool = abs(1 - res) < eps

229 norm_bool *= tmp_bool

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

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

232 return norm_bool

233

234

235@handle_k(mode="skip")

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

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

238

239 Args:

240 obj: Atoms or SCF object.

241 func: A discretized set of functions.

242

243 Keyword Args:

244 eps: Tolerance for the condition.

245

246 Returns:

247 Orthonormality status for the set of functions.

248 """

249 atoms = obj._atoms

250 ortho_bool = check_ortho(atoms, func, eps)

251 norm_bool = check_norm(atoms, func, eps)

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

253 return ortho_bool * norm_bool

254

255

256def get_isovalue(n, percent=85):

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

258

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

260

261 Args:

262 n: Real-space electronic density.

263

264 Keyword Args:

265 percent: Amount of density that should be contained.

266

267 Returns:

268 Isovalue that contains the specified percentage of the density.

269 """

270

271 def deviation(isovalue):

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

273 n_mask = np.sum(n[n > isovalue])

274 return abs(percent - (n_mask / n_ref) * 100)

275

276 # Integrated density

277 n_ref = np.sum(n)

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

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

280 res = minimize_scalar(deviation, bounds=(0, np.max(n)), method="bounded")

281 return res.x

282

283

284def get_tautf(scf):

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

286

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

288

289 Args:

290 scf: SCF object.

291

292 Returns:

293 Real-space Thomas-Fermi kinetic energy density.

294 """

295 atoms = scf.atoms

296 # Use the definition with a division by two

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

298

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

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

301 return tautf

302

303

304def get_tauw(scf):

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

306

307 Reference: Z. Phys. 96, 431.

308

309 Args:

310 scf: SCF object.

311

312 Returns:

313 Real-space von Weizsaecker kinetic energy density.

314 """

315 atoms = scf.atoms

316 if scf.dn_spin is None:

317 dn_spin = get_grad_field(atoms, scf.n_spin)

318 else:

319 dn_spin = scf.dn_spin

320 dn2 = norm(dn_spin, axis=2) ** 2

321 # Use the definition with a division by two

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

323

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

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

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

327 return tauw

328

329

330def get_elf(scf):

331 """Calculate the electron localization function.

332

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

334

335 Args:

336 scf: SCF object.

337

338 Returns:

339 Real-space electron localization function.

340 """

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

342 D0 = get_tautf(scf)

343 X = D / D0

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

345

346

347def get_reduced_gradient(scf, eps=0):

348 """Calculate the reduced density gradient s.

349

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

351

352 Args:

353 scf: SCF object.

354

355 Kwargs:

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

357

358 Returns:

359 Real-space reduced density gradient.

360 """

361 atoms = scf.atoms

362 if scf.dn_spin is None:

363 dn_spin = get_grad_field(atoms, scf.n_spin)

364 else:

365 dn_spin = scf.dn_spin

366 norm_dn = norm(np.sum(dn_spin, axis=0), axis=1)

367

368 kf = (3 * np.pi**2 * scf.n) ** (1 / 3)

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

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

371 s[scf.n < eps] = 0

372 return s

373

374

375def get_spin_squared(scf):

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

377

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

379

380 Args:

381 scf: SCF object.

382

383 Returns:

384 The DFT value for <S^2>.

385 """

386 atoms = scf.atoms

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

388 if not atoms.unrestricted:

389 return 0

390

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

392 rhoXr[rhoXr < 0] = 0

393 rhoX = np.sum(rhoXr) * atoms.dV

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

395 return SX * (SX + 1) + rhoX

396

397

398def get_multiplicity(scf):

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

400

401 Args:

402 scf: SCF object.

403

404 Returns:

405 Multiplicity 2S+1.

406 """

407 S2 = get_spin_squared(scf)

408 # <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

409 S = np.sqrt(S2 + 0.25) - 0.5

410 return 2 * S + 1

411

412

413def get_magnetization(scf):

414 """Calculate the total magnetization M.

415

416 Args:

417 scf: SCF object.

418

419 Returns:

420 Total magnetization.

421 """

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

423 if not scf.atoms.unrestricted:

424 return 0

425

426 return np.sum(scf.n_spin[0] - scf.n_spin[1]) / np.sum(scf.n)

427

428

429def get_bandgap(scf):

430 """Calculate the band gap.

431

432 Args:

433 scf: SCF object.

434

435 Returns:

436 Band gap energy.

437 """

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

439

440 if scf.Z is None:

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

442 return 0

443

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

445 return np.min(e_unocc) - np.max(e_occ)

446

447

448def get_Efermi(obj, epsilon=None):

449 """Calculate the Fermi energy.

450

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

452

453 Args:

454 obj: SCF or Occupations object.

455

456 Keyword Args:

457 epsilon: Eigenenergies.

458

459 Returns:

460 Fermi energy.

461 """

462 # Handle the obj argument

463 if hasattr(obj, "smearing"):

464 if epsilon is None:

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

466 occ = obj

467 else:

468 occ = obj.atoms.occ

469

470 # Calculate the eigenenergies if necessary

471 if epsilon is None:

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

473 else:

474 e_occ = epsilon

475

476 def electron_root(Efermi):

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

478 occ_sum = 0

479 for ik in range(occ.Nk):

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

481 return occ_sum * 2 / occ.Nspin - occ.Nelec

482

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

484 if occ.smearing > 0:

485 return root_scalar(electron_root, bracket=(np.min(e_occ), np.max(e_occ))).root

486

487 if obj.Z is None:

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

489 return np.max(e_occ)

490

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

492 return np.max(e_occ) + (np.min(e_unocc) - np.max(e_occ)) / 2

493

494

495def fermi_distribution(E, mu, kbT):

496 """Calculate the Fermi distribution.

497

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

499

500 Args:

501 E: State energy.

502 mu: Chemical energy or Fermi energy.

503 kbT: Thermic energy or smearing width.

504

505 Returns:

506 Fermi distribution.

507 """

508 x = (E - mu) / kbT

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

510 return 1 / (np.exp(x) + 1)

511

512

513def electronic_entropy(E, mu, kbT):

514 """Calculate the electronic entropic energy.

515

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

517

518 Args:

519 E: State energy.

520 mu: Chemical energy or Fermi energy.

521 kbT: Thermic energy or smearing width.

522

523 Returns:

524 Electronic entropic energy.

525 """

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

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

528 return 0

529 f = fermi_distribution(E, mu, kbT)

530 return f * np.log(f) + (1 - f) * np.log(1 - f)

531

532

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

534 """Calculate the total density of states.

535

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

537

538 Args:

539 epsilon: Eigenenergies.

540 wk: Chemical energy or Fermi energy.

541

542 Keyword Args:

543 spin: Spin channel.

544 npts: Number of energy discretizations.

545 width: Gaussian width.

546

547 Returns:

548 Eigenenergies and DOS.

549 """

550

551 def delta(x, x0, width):

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

553 return np.exp(-(((x - x0) / width) ** 2)) / (np.sqrt(np.pi) * width)

554

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

556 emin = np.min(energies) - 5 * width

557 emax = np.max(energies) + 5 * width

558 e = np.linspace(emin, emax, npts)

559 dos_e = np.zeros(npts)

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

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

562 return e, dos_e