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

1r"""This module calculates the orbital magnetic moment vector for each atom. 

2 

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): 

6 

7:: 

8 

9 === 

10 a e \ a a 

11 m = - -- / D L 

12 orb,v 2m === ii' vii' 

13 ii' 

14 

15with L^a_vii' containing the matrix elements of the angular momentum operator 

16between two partial waves centred at atom a. 

17 

18The orbital magnetic moments are returned in units of μ_B without the sign of 

19the negative electronic charge, q = - e. 

20""" 

21 

22import numpy as np 

23from gpaw.new import zips 

24from gpaw.spinorbit import get_L_vlmm 

25 

26L_vlmm = get_L_vlmm() 

27 

28 

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. 

32 

33 This method assumes that D_asii is on every rank and not parallelised. 

34 

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 

44 

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. 

48 

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. 

53 

54 Partials with unbounded radial functions (negative n_j) are skipped. 

55 """ 

56 

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 

61 

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