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