Coverage for gpaw/new/orbmag.py: 100%
19 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-08 00:17 +0000
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-08 00:17 +0000
1r"""This module calculates the orbital magnetic moment vector for each atom.
3The orbital magnetic moment is calculated in the atom-centred approximation
4where only the PAW correction to the wave function is assumed to contribute.
5This leads to the equation (presented in SI units):
7::
9 ===
10 a e \ a a
11 m = - -- / D L
12 orb,v 2m === ii' vii'
13 ii'
15with L^a_vii' containing the matrix elements of the angular momentum operator
16between two partial waves centred at atom a.
18The orbital magnetic moments are returned in units of μ_B without the sign of
19the negative electronic charge, q = - e.
20"""
22import numpy as np
23from gpaw.new import zips
24from gpaw.spinorbit import get_L_vlmm
26L_vlmm = get_L_vlmm()
29def calculate_orbmag_from_density(D_asii, n_aj, l_aj):
30 """Returns orbital magnetic moment vectors for each atom a
31 calculated from its respective atomic density matrix.
33 This method assumes that D_asii is on every rank and not parallelised.
35 Parameters
36 ----------
37 D_asii : AtomArrays or dictionary
38 Atomic density matrix for each atom a. The i-index refers to the
39 partial waves of an atom and the s-index refers to 0, x, y, and z.
40 n_aj : List of lists of integers
41 Principal quantum number for each radial partial wave j
42 l_aj : List of lists of integers
43 Angular momentum quantum number for each radial partial wave j
45 NB: i is an index for all partial waves for one atom and j is an index for
46 only the radial wave function which is used to build all of the partial
47 waves. i and j do not refer to the same kind of index.
49 Only pairs of partial waves with the same radial function may yield
50 nonzero contributions. The sum can therefore be limited to diagonal blocks
51 of shape [2 * l_j + 1, 2 * l_j + 1] where l_j is the angular momentum
52 quantum number of the j'th radial function.
54 Partials with unbounded radial functions (negative n_j) are skipped.
55 """
57 orbmag_av = np.zeros([len(n_aj), 3])
58 for (a, D_sii), n_j, l_j in zips(D_asii.items(), n_aj, l_aj):
59 assert D_sii.shape[0] == 4
60 D_ii = D_sii[0] # Only the electron density
62 Ni = 0
63 for n, l in zips(n_j, l_j):
64 Nm = 2 * l + 1
65 if n < 0:
66 Ni += Nm
67 continue
68 for v in range(3):
69 orbmag_av[a, v] += np.einsum('ij,ij->',
70 D_ii[Ni:Ni + Nm, Ni:Ni + Nm],
71 L_vlmm[v][l]).real
72 Ni += Nm
73 return orbmag_av