Coverage for gpaw/hubbard.py: 86%
115 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1from typing import List, Tuple
3import numpy as np
4from ase.units import Ha
6from gpaw.typing import Array1D, Array2D, Array3D, ArrayLike2D
9def parse_hubbard_string(type: str) -> Tuple[str, 'HubbardU']:
11 # Parse DFT+U parameters from type-string:
12 # Examples: "type:l,U" or "type:l,U,scale"
13 type, lus = type.split(':')
14 if type == '':
15 type = 'paw'
17 l = []
18 U = []
19 scale = []
21 for lu in lus.split(';'): # Multiple U corrections
22 l_, u_, scale_ = (lu + ',,').split(',')[:3]
23 l.append('spdf'.find(l_))
24 U.append(float(u_) / Ha)
25 if scale_:
26 scale.append(bool(int(scale_)))
27 else:
28 scale.append(True)
29 return type, HubbardU(U, l, scale)
32class HubbardU:
33 def __init__(self, U, l, scale=1):
34 self.scale = scale
35 self.U = U
36 self.l = l
38 def _tuple(self):
39 # Tests use this method to compare to expected values
40 return (self.l, self.U, self.scale)
42 def calculate(self, setup, D_sii):
43 eU = 0.
44 dHU_sii = np.zeros_like(D_sii)
45 for l, U, scale in zip(self.l, self.U, self.scale):
46 nl = np.argwhere(np.equal(setup.l_j, l))
47 if not (len(nl) == 1 or len(nl) == 2):
48 raise RuntimeError(f'Setup has {len(nl)} radial partial waves '
49 f'with angular momentum quantum number {l}.'
50 ' Must be 1 or 2 for DFT+U.')
51 if (len(nl) == len(np.argwhere(np.equal(np.array(setup.n_j)[nl],
52 -1))) and scale == 1):
53 raise RuntimeError('DFT+U correction cannot be scaled if '
54 'there is no bounded partial waves.')
56 eU1, dHU1_sii = hubbard(D_sii, U=U, l=l,
57 l_j=setup.l_j, n_j=setup.n_j,
58 N0_q=setup.N0_q, scale=scale)
59 eU += eU1.real
60 dHU_sii += dHU1_sii
61 return eU, dHU_sii
63 def descriptions(self):
64 for U, l, scale in zip(self.U, self.l, self.scale):
65 yield f'Hubbard: {{U: {U * Ha}, # eV\n'
66 yield f' l: {l},\n'
67 yield f' scale: {bool(scale)}}}'
70def hubbard(D_sii: Array3D,
71 U: float,
72 l: int,
73 l_j: List[int],
74 n_j: List[int],
75 N0_q: Array1D,
76 scale: bool) -> Tuple[float, ArrayLike2D]:
77 nspins = len(D_sii)
79 # j-indices which have the correct angular momentum quantum number
80 nl = np.where(np.equal(l_j, l))[0]
82 nm_j = 2 * np.array(l_j) + 1
83 nm = nm_j[nl[0]]
85 # Get relevant entries of the density matrix
86 i1 = slice(nm_j[:nl[0]].sum(), nm_j[:nl[0]].sum() + nm)
88 eU = 0.
89 dHU_sii = np.zeros_like(D_sii)
91 for s, D_ii in enumerate(D_sii):
92 N_mm, dHU_ii = aoom(D_ii, l, l_j, n_j, N0_q, scale)
93 N_mm = N_mm / 2 * nspins
95 if nspins == 4:
96 N_mm = N_mm / 2.0
97 if s == 0:
98 eU1 = U / 2. * (N_mm - 0.5 * N_mm @ N_mm).trace()
100 dHU_mm = U / 2. * (np.eye(nm) - N_mm.T)
102 else:
103 eU1 = -U / 2. * (0.5 * N_mm @ N_mm).trace()
105 dHU_mm = -U / 2. * N_mm.T
106 else:
107 eU1 = U / 2. * (N_mm - N_mm @ N_mm).trace()
109 dHU_mm = U / 2. * (np.eye(nm) - 2 * N_mm)
111 eU += eU1
112 if nspins == 1:
113 # Add contribution from other spin manifold
114 eU += eU1
116 if len(nl) == 1:
117 dHU_ii[i1, i1] *= dHU_mm
118 elif len(nl) == 2:
119 i2 = slice(nm_j[:nl[1]].sum(), nm_j[:nl[1]].sum() + nm)
121 dHU_ii[i1, i1] *= dHU_mm
122 dHU_ii[i1, i2] *= dHU_mm
123 dHU_ii[i2, i1] *= dHU_mm
124 dHU_ii[i2, i2] *= dHU_mm
125 else:
126 raise NotImplementedError
128 dHU_sii[s] = dHU_ii
130 return eU, dHU_sii
133def aoom(D_ii: Array2D,
134 l: int,
135 l_j: List[int],
136 n_j: List[int],
137 N0_q: Array1D,
138 scale: bool = True) -> Tuple[Array2D, Array2D]:
139 """
140 This function returns the atomic orbital occupation matrix (aoom) for a
141 given l quantum number.
143 The submatrix / submatrices of the density matrix (D_ii) for the
144 selected l quantum number are determined and summed together which
145 represents the orbital occupation matrix. For l=2 this is a 5x5 matrix.
147 If scale = True, the inner products are scaled such that the inner product
148 of the bounded partial waves is 1.
149 """
151 # j-indices which have the correct angular momentum quantum number
152 nl = np.where(np.equal(l_j, l))[0]
154 nm_j = 2 * np.array(l_j) + 1
155 nm = nm_j[nl[0]]
157 # Get relevant entries of the density matrix
158 i1 = slice(nm_j[:nl[0]].sum(), nm_j[:nl[0]].sum() + nm)
160 dHU_ii = np.zeros_like(D_ii)
161 if len(nl) == 2:
162 # First get q-indices for the radial inner products
163 q1 = nl[0] * len(l_j) - (nl[0] - 1) * nl[0] // 2 # Bounded-bounded
164 q2 = nl[1] * len(l_j) - (nl[1] - 1) * nl[1] // 2 # Unbounded-unbounded
165 q12 = q1 + nl[1] - nl[0] # Bounded-unbounded
167 # If the Hubbard correction should be scaled, the three inner products
168 # will be divided by the inner product of the bounded partial wave,
169 # increasing the value of these inner products since 0 < N0_q[q1] < 1
170 if scale:
171 if n_j[nl[1]] == -1:
172 N1 = 1
173 N2 = N0_q[q2] / N0_q[q1]
174 N12 = N0_q[q12] / N0_q[q1]
175 else:
176 N1 = 1
177 N2 = 1
178 N12 = N0_q[q12] / np.sqrt(N0_q[q1] * N0_q[q2])
179 else:
180 N1 = N0_q[q1]
181 N2 = N0_q[q2]
182 N12 = N0_q[q12]
184 # Get the entries in the density matrix of the unbounded partial waves
185 i2 = slice(nm_j[:nl[1]].sum(), nm_j[:nl[1]].sum() + nm)
187 # Scale and add the four submatrices to get the occupation matrix
188 N_mm = D_ii[i1, i1] * N1 + D_ii[i2, i2] * N2 + (D_ii[i1, i2]
189 + D_ii[i2, i1]) * N12
191 dHU_ii[i1, i1] = N1
192 dHU_ii[i1, i2] = N12
193 dHU_ii[i2, i1] = N12
194 dHU_ii[i2, i2] = N2
196 return N_mm, dHU_ii
197 elif len(nl) == 1:
198 q1 = nl[0] * len(l_j) - (nl[0] - 1) * nl[0] // 2
199 if scale:
200 N1 = 1
201 else:
202 N1 = N0_q[q1]
204 N_mm = D_ii[i1, i1] * N1
205 dHU_ii[i1, i1] = N1
206 return N_mm, dHU_ii
207 else:
208 raise NotImplementedError(f'Setup has {len(nl)} partial waves with '
209 f'angular momentum quantum number {l}. '
210 'Must be 1 or 2 for DFT+U correction.')