Coverage for eminus/tools.py: 97.46%

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 masses is None:

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

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

60@handle_k

61def orbital_center(obj, psirs):

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

64 Args:

65 obj: Atoms or SCF object.

66 psirs: Set of orbitals in real-space.

68 Returns:

69 Center of masses.

70 """

71 atoms = obj._atoms

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

74 Ncom = psirs.shape[2]

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

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

78 # Squared orbitals

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

80 for i in range(Ncom):

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

82 coms[spin] = coms_spin

83 return coms

86def inertia_tensor(coords, masses=None):

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

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

91 Args:

92 coords: Array of real-space coordinates.

94 Keyword Args:

95 masses: Mass or weight for each coordinate.

97 Returns:

98 Inertia tensor.

99 """

100 if masses is None:

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

102

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

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

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

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

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

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

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

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

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

112 return I

113

114

115def get_dipole(scf, n=None):

116 """Calculate the electric dipole moment.

117

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

119

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

121

122 Args:

123 scf: SCF object.

124

125 Keyword Args:

126 n: Real-space electronic density.

127

128 Returns:

129 Electric dipole moment in e * Bohr.

130 """

131 atoms = scf.atoms

132 if n is None:

133 if scf.n is None:

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

135 return 0

136 n = scf.n

137

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

139 mu = np.zeros(3)

140 for i in range(atoms.Natoms):

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

142

143 for dim in range(3):

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

145 return mu

146

147

148def get_ip(scf):

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

150

151 Reference: Physica 1, 104.

152

153 Args:

154 scf: SCF object.

155

156 Returns:

157 Ionization potential in Hartree.

158 """

159 scf.kpts._assert_gamma_only()

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

161 # Account for spin-polarized calculations

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

163 return -epsilon[-1]

164

165

166@handle_k(mode="skip")

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

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

169

170 Args:

171 obj: Atoms or SCF object.

172 func: A discretized set of functions.

173

174 Keyword Args:

175 eps: Tolerance for the condition.

176

177 Returns:

178 Orthogonality status for the set of functions.

179 """

180 atoms = obj._atoms

181 func = np.atleast_3d(func)

182

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

184 if atoms.occ.Nstate == 1:

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

186 return True

187

188 ortho_bool = True

189 # Check the condition for every combination

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

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

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

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

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

195 tmp_bool = abs(res) < eps

196 ortho_bool *= tmp_bool

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

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

199 return ortho_bool

200

201

202@handle_k(mode="skip")

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

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

205

206 Args:

207 obj: Atoms or SCF object.

208 func: A discretized set of functions.

209

210 Keyword Args:

211 eps: Tolerance for the condition.

212

213 Returns:

214 Normalization status for the set of functions.

215 """

216 atoms = obj._atoms

217 func = np.atleast_3d(func)

218

219 norm_bool = True

220 # Check the condition for every function

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

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

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

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

225 tmp_bool = abs(1 - res) < eps

226 norm_bool *= tmp_bool

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

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

229 return norm_bool

230

231

232@handle_k(mode="skip")

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

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

235

236 Args:

237 obj: Atoms or SCF object.

238 func: A discretized set of functions.

239

240 Keyword Args:

241 eps: Tolerance for the condition.

242

243 Returns:

244 Orthonormality status for the set of functions.

245 """

246 atoms = obj._atoms

247 ortho_bool = check_ortho(atoms, func, eps)

248 norm_bool = check_norm(atoms, func, eps)

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

250 return ortho_bool * norm_bool

251

252

253def get_isovalue(n, percent=85):

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

255

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

257

258 Args:

259 n: Real-space electronic density.

260

261 Keyword Args:

262 percent: Amount of density that should be contained.

263

264 Returns:

265 Isovalue that contains the specified percentage of the density.

266 """

267

268 def deviation(isovalue):

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

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

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

272

273 # Integrated density

274 n_ref = np.sum(n)

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

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

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

278 return res.x

279

280

281def get_tautf(scf):

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

283

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

285

286 Args:

287 scf: SCF object.

288

289 Returns:

290 Real-space Thomas-Fermi kinetic energy density.

291 """

292 atoms = scf.atoms

293 # Use the definition with a division by two

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

295

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

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

298 return tautf

299

300

301def get_tauw(scf):

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

303

304 Reference: Z. Phys. 96, 431.

305

306 Args:

307 scf: SCF object.

308

309 Returns:

310 Real-space von Weizsaecker kinetic energy density.

311 """

312 atoms = scf.atoms

313 if scf.dn_spin is None:

314 dn_spin = get_grad_field(atoms, scf.n_spin)

315 else:

316 dn_spin = scf.dn_spin

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

318 # Use the definition with a division by two

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

320

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

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

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

324 return tauw

325

326

327def get_elf(scf):

328 """Calculate the electron localization function.

329

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

331

332 Args:

333 scf: SCF object.

334

335 Returns:

336 Real-space electron localization function.

337 """

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

339 D0 = get_tautf(scf)

340 X = D / D0

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

342

343

344def get_reduced_gradient(scf, eps=0):

345 """Calculate the reduced density gradient s.

346

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

348

349 Args:

350 scf: SCF object.

351

352 Kwargs:

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

354

355 Returns:

356 Real-space reduced density gradient.

357 """

358 atoms = scf.atoms

359 if scf.dn_spin is None:

360 dn_spin = get_grad_field(atoms, scf.n_spin)

361 else:

362 dn_spin = scf.dn_spin

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

364

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

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

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

368 s[scf.n < eps] = 0

369 return s

370

371

372def get_spin_squared(scf):

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

374

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

376

377 Args:

378 scf: SCF object.

379

380 Returns:

381 The DFT value for <S^2>.

382 """

383 atoms = scf.atoms

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

385 if not atoms.unrestricted:

386 return 0

387

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

389 rhoXr[rhoXr < 0] = 0

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

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

392 return SX * (SX + 1) + rhoX

393

394

395def get_multiplicity(scf):

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

397

398 Args:

399 scf: SCF object.

400

401 Returns:

402 Multiplicity 2S+1.

403 """

404 S2 = get_spin_squared(scf)

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

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

407 return 2 * S + 1

408

409

410def get_magnetization(scf):

411 """Calculate the total magnetization M.

412

413 Args:

414 scf: SCF object.

415

416 Returns:

417 Total magnetization.

418 """

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

420 if not scf.atoms.unrestricted:

421 return 0

422

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

424

425

426def get_bandgap(scf):

427 """Calculate the band gap.

428

429 Args:

430 scf: SCF object.

431

432 Returns:

433 Band gap energy.

434 """

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

436

437 if scf.Z is None:

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

439 return 0

440

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

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

443

444

445def get_Efermi(obj, epsilon=None):

446 """Calculate the Fermi energy.

447

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

449

450 Args:

451 obj: SCF or Occupations object.

452

453 Keyword Args:

454 epsilon: Eigenenergies.

455

456 Returns:

457 Fermi energy.

458 """

459 # Handle the obj argument

460 if hasattr(obj, "smearing"):

461 if epsilon is None:

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

463 occ = obj

464 else:

465 occ = obj.atoms.occ

466

467 # Calculate the eigenenergies if necessary

468 if epsilon is None:

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

470 else:

471 e_occ = epsilon

472

473 def electron_root(Efermi):

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

475 occ_sum = 0

476 for ik in range(occ.Nk):

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

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

479

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

481 if occ.smearing > 0:

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

483

484 if obj.Z is None:

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

486 return np.max(e_occ)

487

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

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

490

491

492def fermi_distribution(E, mu, kbT):

493 """Calculate the Fermi distribution.

494

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

496

497 Args:

498 E: State energy.

499 mu: Chemical energy or Fermi energy.

500 kbT: Thermic energy or smearing width.

501

502 Returns:

503 Fermi distribution.

504 """

505 x = (E - mu) / kbT

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

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

508

509

510def electronic_entropy(E, mu, kbT):

511 """Calculate the electronic entropic energy.

512

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

514

515 Args:

516 E: State energy.

517 mu: Chemical energy or Fermi energy.

518 kbT: Thermic energy or smearing width.

519

520 Returns:

521 Electronic entropic energy.

522 """

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

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

525 return 0

526 f = fermi_distribution(E, mu, kbT)

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

528

529

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

531 """Calculate the total density of states.

532

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

534

535 Args:

536 epsilon: Eigenenergies.

537 wk: Chemical energy or Fermi energy.

538

539 Keyword Args:

540 spin: Spin channel.

541 npts: Number of energy discretizations.

542 width: Gaussian width.

543

544 Returns:

545 Eigenenergies and DOS.

546 """

547

548 def delta(x, x0, width):

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

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

551

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

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

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

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

556 dos_e = np.zeros(npts)

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

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

559 return e, dos_e