Coverage for tbkit/meanfield.py: 100%
127 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-30 11:04 +0000
1r"""
2The Hubbard model in the (collinear, unrestricted Hartree-Fock) mean-field
3approximation:
5.. math::
7 H = \sum_{ij\sigma} t_{ij}c^\dagger_{i\sigma}c_{j\sigma}
8 + U\sum_i n_{i\uparrow}n_{i\downarrow}
9 \;\to\;
10 \sum_\sigma\Big[\sum_{ij} t_{ij}c^\dagger_{i\sigma}c_{j\sigma}
11 + U\sum_i\langle n_{i\bar\sigma}\rangle n_{i\sigma}\Big]
12 - U\sum_i\langle n_{i\uparrow}\rangle\langle n_{i\downarrow}\rangle\, .
14Each spin moves in the density of the other; the densities are iterated
15(with linear mixing) until self-consistent. With :math:`U > 0` on a
16bipartite lattice at half filling it reproduces Lieb's theorem (the total
17moment :math:`S = ||A| - |B||/2`) and the edge magnetism of zigzag graphene.
18"""
19from __future__ import annotations
21import numpy as np
22from numpy.typing import ArrayLike, NDArray
23import scipy.linalg as LA
25import tbkit.error_handling as error_handling
26import tbkit.occupation as occupation
29class MeanFieldResult():
30 '''
31 Self-consistent solution of *hubbard_mean_field*.
33 Attributes
34 ----------
35 n_up : Real ndarray
36 Spin-up density on every site.
37 n_dn : Real ndarray
38 Spin-down density on every site.
39 magnetization : Real ndarray
40 :math:`m_i = (n_{i\\uparrow} - n_{i\\downarrow})/2`.
41 en_up : Real ndarray
42 Spin-up mean-field energies.
43 en_dn : Real ndarray
44 Spin-down mean-field energies.
45 e_fermi : Real number
46 Fermi level.
47 energy : Real number
48 Mean-field ground-state (free, at T > 0) energy
49 :math:`\\sum f E - U\\sum_i n_{i\\uparrow}n_{i\\downarrow}`.
50 iterations : Integer
51 Number of iterations.
52 '''
54 def __init__(self, n_up, n_dn, en_up, en_dn, e_fermi, energy, iterations) -> None:
55 self.n_up, self.n_dn = n_up, n_dn
56 self.magnetization = (n_up - n_dn) / 2
57 self.en_up, self.en_dn = en_up, en_dn
58 self.e_fermi = e_fermi
59 self.energy = energy
60 self.iterations = iterations
62 @property
63 def total_spin(self) -> float:
64 r'''
65 Total :math:`S_z = \sum_i m_i`.
66 '''
67 return float(np.sum(self.magnetization))
70def _occupy(en_up, en_dn, n_electrons, temperature):
71 '''
72 Private function. Occupations of both spins sharing one Fermi level: at
73 T = 0 the lowest levels, the electrons left for a degenerate level at
74 the Fermi energy shared equally among its states (a fixed, arbitrary
75 choice among them would make the iteration oscillate), else Fermi-Dirac.
76 '''
77 energies = np.concatenate([en_up, en_dn])
78 if temperature == 0:
79 n = int(round(n_electrons))
80 order = np.sort(energies)
81 f = np.zeros(len(energies))
82 if n:
83 e_f = order[n - 1]
84 tol = 1e-8 * max(1., np.max(np.abs(energies)))
85 below = energies < e_f - tol
86 shell = np.abs(energies - e_f) <= tol
87 f[below] = 1.
88 f[shell] = (n - below.sum()) / shell.sum()
89 mu = occupation.fermi_level(energies, n_electrons)
90 else:
91 mu = occupation.fermi_level(energies, n_electrons, temperature)
92 f = occupation.fermi_dirac(energies, mu, temperature)
93 return f[:len(en_up)], f[len(en_up):], mu
96def hubbard_mean_field(
97 ham, U: float | ArrayLike, n_electrons: float, temperature: float = 0.,
98 n_up: ArrayLike | None = None, n_dn: ArrayLike | None = None, mixing: float = 0.5,
99 tol: float = 1e-8, max_iter: int = 20000, seed=None,
100) -> MeanFieldResult:
101 r'''
102 Solve the Hubbard model in mean field, self-consistently.
104 The iteration converges to *a* self-consistent solution, not
105 necessarily the one of lowest energy: compare the *energy* of the
106 solutions reached from several starts (random seeds, or a staggered
107 magnetization on the two sublattices of a bipartite lattice).
109 Parameters
110 ----------
111 ham
112 Square Hermitian matrix (sparse or dense), the spin-independent
113 single-particle Hamiltonian (e.g. *System.ham*).
114 U : Real number, or array of one per site
115 Onsite repulsion.
116 n_electrons : Positive real, at most 2N
117 Number of electrons (both
118 spins; N for half filling). An integer at T = 0.
119 temperature : Positive real or zero, default 0
120 n_up : Real array, default None
121 Initial spin-up density.
122 n_dn : Real array, default None
123 Initial spin-down density.
124 By default, a small random magnetization on top of the uniform
125 density (to break the spin symmetry).
126 mixing : Real in (0, 1], default 0.5
127 Fraction of the new
128 density mixed in at each iteration. (Plain linear mixing converges
129 only to stable solutions, but slowly near a magnetic transition:
130 raise *max_iter* there.)
131 tol : Positive real, default 1e-8
132 Convergence threshold on
133 the largest density change.
134 max_iter : Positive integer, default 20000
135 seed
136 Default value None. Seed of the random initial magnetization.
138 Returns
139 -------
140 result : :class:`MeanFieldResult`
141 '''
142 ham = ham.toarray() if hasattr(ham, 'toarray') else np.asarray(ham)
143 error_handling.square_matrix(ham, 'ham')
144 n = len(ham)
145 error_handling.hermitian_dense(ham)
146 U = np.broadcast_to(np.asarray(U, dtype='f8'), (n,)).copy() if np.ndim(U) == 0 \
147 else np.asarray(U, dtype='f8')
148 error_handling.ndarray(U, 'U', n)
149 error_handling.electrons(n_electrons, 2 * n)
150 error_handling.positive_real_zero(temperature, 'temperature')
151 if temperature == 0:
152 error_handling.integer_electrons(n_electrons)
153 error_handling.mixing(mixing)
154 error_handling.positive_real(tol, 'tol')
155 error_handling.positive_int(max_iter, 'max_iter')
156 if n_up is None or n_dn is None:
157 rng = np.random.default_rng(seed)
158 kick = 0.1 * rng.uniform(-1., 1., n)
159 n_up = np.full(n, n_electrons / (2 * n)) + kick
160 n_dn = np.full(n, n_electrons / (2 * n)) - kick
161 up, dn = np.asarray(n_up, dtype='f8').copy(), np.asarray(n_dn, dtype='f8').copy()
162 error_handling.ndarray(up, 'n_up', n)
163 error_handling.ndarray(dn, 'n_dn', n)
164 for it in range(1, max_iter + 1):
165 en_up, v_up = LA.eigh(ham + np.diag(U * dn))
166 en_dn, v_dn = LA.eigh(ham + np.diag(U * up))
167 f_up, f_dn, mu = _occupy(en_up, en_dn, n_electrons, temperature)
168 new_up = (np.abs(v_up) ** 2) @ f_up
169 new_dn = (np.abs(v_dn) ** 2) @ f_dn
170 change = max(np.max(np.abs(new_up - up)), np.max(np.abs(new_dn - dn)))
171 up = (1 - mixing) * up + mixing * new_up
172 dn = (1 - mixing) * dn + mixing * new_dn
173 if change < tol:
174 break
175 else:
176 error_handling.converged(False, 'hubbard_mean_field')
177 energy = float(np.sum(f_up * en_up) + np.sum(f_dn * en_dn) - np.sum(U * new_up * new_dn))
178 return MeanFieldResult(new_up, new_dn, en_up, en_dn, mu, energy, it)
181#: Pauli matrices sigma_x, sigma_y, sigma_z.
182_SIGMA = np.array([[[0., 1.], [1., 0.]], [[0., -1j], [1j, 0.]], [[1., 0.], [0., -1.]]], dtype='c16')
185class NonCollinearResult():
186 r'''
187 Self-consistent solution of *hubbard_mean_field_noncollinear*.
189 Attributes
190 ----------
191 rho : Complex ndarray, shape (N, 2, 2)
192 Onsite spin density
193 matrices :math:`\rho_{i,\alpha\beta} = \langle c^\dagger_{i\beta}c_{i\alpha}\rangle`.
194 density : Real ndarray
195 :math:`n_i = \mathrm{Tr}\,\rho_i`.
196 magnetization : Real ndarray, shape (N, 3)
197 :math:`\mathbf{m}_i = \mathrm{Tr}(\rho_i\boldsymbol\sigma)/2`.
198 energies : Real ndarray
199 Mean-field energies (2N).
200 states : Complex ndarray, shape (2N, 2N)
201 Eigenvectors (columns),
202 rows site-major (spin up, down).
203 e_fermi : Real number
204 Fermi level.
205 energy : Real number
206 Mean-field ground-state energy
207 :math:`\sum f E - U\sum_i(\rho_{i\uparrow\uparrow}\rho_{i\downarrow\downarrow}
208 - |\rho_{i\uparrow\downarrow}|^2)`.
209 iterations : Integer
210 Number of iterations.
211 '''
213 def __init__(self, rho, energies, states, e_fermi, energy, iterations) -> None:
214 self.rho = rho
215 self.density = np.real(np.trace(rho, axis1=1, axis2=2))
216 self.magnetization = np.real(np.einsum('iab,sba->is', rho, _SIGMA)) / 2
217 self.energies, self.states = energies, states
218 self.e_fermi = e_fermi
219 self.energy = energy
220 self.iterations = iterations
222 @property
223 def total_spin(self) -> NDArray[np.float64]:
224 r'''
225 Total spin :math:`\sum_i\mathbf{m}_i`, shape (3,).
226 '''
227 return self.magnetization.sum(axis=0)
230def hubbard_mean_field_noncollinear(
231 ham, U: float | ArrayLike, n_electrons: float, temperature: float = 0.,
232 magnetization: ArrayLike | None = None, spinful: bool = False, mixing: float = 0.5,
233 tol: float = 1e-8, max_iter: int = 20000, seed=None,
234) -> NonCollinearResult:
235 r'''
236 Solve the Hubbard model in the non-collinear (spin-rotation invariant)
237 Hartree-Fock approximation: the onsite spin density matrix
238 :math:`\rho_i` is kept whole, so the moments may point in any
239 direction (spirals, the 120-degree order of frustrated lattices,
240 moments canted by spin-orbit coupling). The Hartree and Fock terms of
241 :math:`Un_{i\uparrow}n_{i\downarrow}` add, on every site,
243 .. math::
245 V_i = U\,(n_i\mathbb{1} - \rho_i) = U\begin{pmatrix}\rho_{\downarrow\downarrow} & -\rho_{\uparrow\downarrow}\\
246 -\rho_{\downarrow\uparrow} & \rho_{\uparrow\uparrow}\end{pmatrix}
247 = U\Big(\frac{n_i}{2}\mathbb{1} - \mathbf{m}_i\cdot\boldsymbol\sigma\Big)\, ,
249 iterated with linear mixing until :math:`\rho` is self-consistent. With
250 :math:`\rho_{\uparrow\downarrow} = 0` it is *hubbard_mean_field*
251 (collinear moments along :math:`z` stay collinear, and give the same
252 solution). As there, several starts should be compared by *energy*.
254 Parameters
255 ----------
256 ham
257 Square Hermitian matrix: the single-particle Hamiltonian,
258 (N, N) spin independent, or, with ``spinful=True``, (2N, 2N) with
259 rows site-major (spin up, down) -- e.g. with spin-orbit coupling.
260 U : Real number, or array of one per site
261 Onsite repulsion.
262 n_electrons : Positive real, at most 2N
263 An integer at T = 0.
264 temperature : Positive real or zero, default 0
265 magnetization : Real array, shape (N, 3)
266 Default value None
267 (random directions of length 0.1). Initial moments, on top of a
268 uniform density.
269 spinful : Boolean, default False
270 See *ham*.
271 mixing : Real in (0, 1], default 0.5
272 tol : Positive real, default 1e-8
273 Convergence threshold on
274 the largest change of :math:`\rho`.
275 max_iter : Positive integer, default 20000
276 seed
277 Default value None. Seed of the random initial moments.
279 Returns
280 -------
281 result : :class:`NonCollinearResult`
282 '''
283 ham = ham.toarray() if hasattr(ham, 'toarray') else np.asarray(ham)
284 error_handling.square_matrix(ham, 'ham')
285 error_handling.hermitian_dense(ham)
286 error_handling.boolean(spinful, 'spinful')
287 if spinful:
288 error_handling.even_dimension(len(ham))
289 h0 = ham.astype('c16')
290 else:
291 h0 = np.kron(ham, np.eye(2)).astype('c16')
292 n = len(h0) // 2
293 U = np.full(n, U, dtype='f8') if np.ndim(U) == 0 else np.asarray(U, dtype='f8')
294 error_handling.ndarray(U, 'U', n)
295 error_handling.electrons(n_electrons, 2 * n)
296 error_handling.positive_real_zero(temperature, 'temperature')
297 if temperature == 0:
298 error_handling.integer_electrons(n_electrons)
299 error_handling.mixing(mixing)
300 error_handling.positive_real(tol, 'tol')
301 error_handling.positive_int(max_iter, 'max_iter')
302 if magnetization is None:
303 vec = np.random.default_rng(seed).normal(size=(n, 3))
304 magnetization = 0.1 * vec / np.linalg.norm(vec, axis=1)[:, None]
305 magnetization = np.asarray(magnetization, dtype='f8')
306 error_handling.magnetization(magnetization, n)
307 # rho = n/2 + m.sigma (so that Tr(rho sigma)/2 = m)
308 rho = (n_electrons / (2 * n) * np.eye(2)[None]
309 + np.einsum('is,sab->iab', magnetization, _SIGMA)).astype('c16')
310 blocks = np.arange(n)
311 for it in range(1, max_iter + 1):
312 pot = U[:, None, None] * (np.trace(rho, axis1=1, axis2=2)[:, None, None] * np.eye(2)[None] - rho)
313 h = h0.copy()
314 for a in range(2):
315 for b in range(2):
316 h[2 * blocks + a, 2 * blocks + b] += pot[:, a, b]
317 en, vec = np.linalg.eigh(h)
318 f, _, mu = _occupy(en, en[:0], n_electrons, temperature)
319 full = (vec * f[None, :]) @ vec.conj().T
320 new = full.reshape(n, 2, n, 2)[blocks, :, blocks, :]
321 change = np.max(np.abs(new - rho))
322 rho = (1 - mixing) * rho + mixing * new
323 if change < tol:
324 break
325 else:
326 error_handling.converged(False, 'hubbard_mean_field_noncollinear')
327 double = U * np.real(new[:, 0, 0] * new[:, 1, 1] - np.abs(new[:, 0, 1]) ** 2)
328 energy = float(np.sum(f * en) - np.sum(double))
329 return NonCollinearResult(new, en, vec, mu, energy, it)