Rustb/response/linear/mod.rs
1//! # Linear response: Berry curvature and quantum metric
2//!
3//! ## Theory
4//!
5//! The quantum geometric tensor is
6//!
7//! $$G^{ab}_n(\mathbf{k}) = \sum_{m=\not n}
8//! \frac{\langle\partial_a u_n|u_m\rangle\langle u_m|\partial_b u_n\rangle}
9//! {(E_n-E_m)^2} = g^{ab}_n - \frac{i}{2}\Omega^{ab}_n$$
10//!
11//! In terms of velocity matrix elements $v^\alpha_{nm} = \langle u_n|\partial_\alpha H|u_m\rangle$:
12//!
13//! $$G^{ab}_n(\mathbf{k}) = \sum_{m=\not n} \frac{v^a_{nm} v^b_{mn}}{(E_n-E_m)^2 + \eta^2}$$
14//!
15//! where $\eta$ is a small regularisation width. The real and imaginary
16//! parts give the quantum metric $g^{ab}_n$ and Berry curvature $\Omega^{ab}_n$:
17//!
18//! $$g^{ab}_n = \operatorname{Re} G^{ab}_n, \qquad
19//! \Omega^{ab}_n = -2\operatorname{Im} G^{ab}_n$$
20//!
21//! The **anomalous Hall conductivity** (AHC) at $T=0$ is
22//!
23//! $$\sigma^{xy}_{\text{AHC}}(\mu) = -\frac{e^2}{\hbar}\sum_n
24//! \int_{\text{BZ}} \Theta(\mu-E_n)\Omega^{xy}_n(\mathbf{k})d\mathbf{k}$$
25//!
26//! ## API
27//!
28//! | Method | Path | Formula |
29//! |--------|------|---------|
30//! | `BerryCurvature::berry_curvature_at` | one k-point | band-resolved $\Omega_n^{ab}$ |
31//! | `Model::hall_conductivity` | direct sum or energy cut | $\sigma_{\text{AHC}}(\mu)$ |
32
33use ndarray::prelude::*;
34use ndarray_linalg::*;
35use rayon::prelude::*;
36
37use crate::Gauge;
38use crate::Model;
39use crate::RMatrixData;
40use crate::error::Result;
41use crate::thermodynamics::Occupation;
42
43use super::config::{Integration, Parameters, mesh_array, parameters_occupation, validate_sorted};
44use super::energy_cut::{integrate_fermi_cut_2d, integrate_fermi_cut_3d};
45use super::kernel::quadrature_occupied_geometry_simplex;
46use super::tracking::{build_tetrahedra_3d, build_triangles_2d, global_band_track};
47use super::types::{SIMPLEX_GAP_TOL, VertexKernel};
48
49/// Hall conductivity evaluated on the requested chemical-potential grid.
50#[derive(Clone, Debug, PartialEq)]
51pub struct HallConductivityResult {
52 /// Chemical potentials copied from the input configuration.
53 pub chemical_potentials: Array1<f64>,
54 /// Hall response at every chemical potential.
55 pub conductivity: Array1<f64>,
56}
57
58impl HallConductivityResult {
59 /// Return the scalar value produced by [`Parameters::at_mu`].
60 pub fn single(&self) -> Option<f64> {
61 (self.conductivity.len() == 1).then(|| self.conductivity[0])
62 }
63}
64
65/// Integrate occupation-weighted Berry curvature and quantum metric over the BZ.
66///
67/// ```text
68/// total_berry(μ) = Σ_n ∫_BZ f(E_n, μ) Ω^{ab}_n(k) dk
69/// total_metric(μ) = Σ_n ∫_BZ f(E_n, μ) g^{ab}_n(k) dk
70/// ```
71///
72/// where `Ω_n = −2 Im G_n`, `g_n = Re G_n`, and
73/// `G_n = Σ_{m≠n} K^{ab}_nm / ((E_n−E_m)² + η²)`.
74///
75/// Inside each simplex, `K_nm` and `E_n` are linearly interpolated,
76/// then the kernel is evaluated at degree‑2 symmetric quadrature points.
77///
78/// Values use fractional-coordinate volume. Divide by `det(lat)` for
79/// Cartesian reciprocal-space volume.
80pub(crate) fn integrate_occupied_geometry(
81 all_pts: &[VertexKernel],
82 k_mesh: &Array1<usize>,
83 eta: f64,
84 chemical_potentials: &Array1<f64>,
85 occupation: Occupation,
86) -> (Array1<f64>, Array1<f64>, usize) {
87 let dim = k_mesh.len();
88 let mut total_g = Array1::<f64>::zeros(chemical_potentials.len());
89 let mut total_o = Array1::<f64>::zeros(chemical_potentials.len());
90 let mut unsafe_count = 0usize;
91
92 match dim {
93 2 => {
94 let (nx, ny) = (k_mesh[0], k_mesh[1]);
95 let inv_nx = 1.0 / nx as f64;
96 let inv_ny = 1.0 / ny as f64;
97 for ix in 0..nx {
98 for iy in 0..ny {
99 let sims = build_triangles_2d(ix, iy, nx, ny, inv_nx, inv_ny, all_pts);
100 for sim in &sims {
101 if sim.diag.min_gap < SIMPLEX_GAP_TOL {
102 unsafe_count += 1;
103 }
104 let (g, o) = quadrature_occupied_geometry_simplex(
105 sim,
106 eta,
107 chemical_potentials,
108 occupation,
109 );
110 total_g += &g;
111 total_o += &o;
112 }
113 }
114 }
115 }
116 3 => {
117 let (nx, ny, nz) = (k_mesh[0], k_mesh[1], k_mesh[2]);
118 let inv_nx = 1.0 / nx as f64;
119 let inv_ny = 1.0 / ny as f64;
120 let inv_nz = 1.0 / nz as f64;
121 for ix in 0..nx {
122 for iy in 0..ny {
123 for iz in 0..nz {
124 let sims = build_tetrahedra_3d(
125 ix, iy, iz, nx, ny, nz, inv_nx, inv_ny, inv_nz, all_pts,
126 );
127 for sim in &sims {
128 if sim.diag.min_gap < SIMPLEX_GAP_TOL {
129 unsafe_count += 1;
130 }
131 let (g, o) = quadrature_occupied_geometry_simplex(
132 sim,
133 eta,
134 chemical_potentials,
135 occupation,
136 );
137 total_g += &g;
138 total_o += &o;
139 }
140 }
141 }
142 }
143 }
144 _ => panic!("linear::integrate: only dim=2,3 supported, got dim={dim}"),
145 }
146
147 (total_g, total_o, unsafe_count)
148}
149
150impl<const SPIN: bool, const DIM: usize, R: RMatrixData> Model<SPIN, DIM, R> {
151 /// Evaluate charge or spin Hall conductivity using the unified
152 /// [`Parameters`] configuration.
153 ///
154 /// Reads `kmesh`, `direction` (rank 2), `mu`, `T`, `eta` and `spin` from
155 /// the parameter set; `omega` and `field_symmetry` are ignored. The
156 /// returned array has one value per chemical potential. Direct
157 /// integration reuses the band-resolved Berry curvature for the entire
158 /// chemical-potential grid; energy-cut integration tracks bands between
159 /// simplex vertices before integrating the occupied region.
160 pub fn hall_conductivity(&self, params: &Parameters<DIM>) -> Result<HallConductivityResult> {
161 params.validate_rank2()?;
162 let spin = params.spin;
163 if !SPIN && let Some(direction) = spin {
164 return Err(crate::TbError::SpinNotAllowed(direction));
165 }
166 match params.integration {
167 Integration::Direct | Integration::EnergyCut => {}
168 Integration::Simplex => {
169 return Err(crate::TbError::InvalidResponseParameter {
170 parameter: "integration",
171 message:
172 "hall_conductivity supports Integration::Direct or EnergyCut, not Simplex"
173 .into(),
174 });
175 }
176 }
177 if params.integration == Integration::EnergyCut {
178 validate_sorted(¶ms.mu, "mu")?;
179 if DIM != 2 && DIM != 3 {
180 return Err(crate::TbError::InvalidDimension {
181 dim: DIM,
182 supported: vec![2, 3],
183 });
184 }
185 }
186
187 let k_mesh = mesh_array(¶ms.kmesh);
188 let determinant = self.lat.det()?;
189 let dir_a = params.direction.row(0).to_owned();
190 let dir_b = params.direction.row(1).to_owned();
191 let occupation = parameters_occupation(params);
192 let conductivity = match params.integration {
193 Integration::Direct => {
194 let kvec: Array2<f64> = crate::kpoints::gen_kmesh(&k_mesh)?;
195 let nk = kvec.nrows();
196 // Validation already happened at the entry point; reuse the
197 // unvalidated kernel instead of re-validating per k-point.
198 let band_data: Vec<Result<_>> = kvec
199 .axis_iter(Axis(0))
200 .into_par_iter()
201 .map(|k| self.berry_curvature_at_impl(&k, params))
202 .collect();
203 let band_data = band_data.into_iter().collect::<Result<Vec<_>>>()?;
204 let values: Vec<f64> = params
205 .mu
206 .par_iter()
207 .map(|&mu| {
208 let sum: f64 = band_data
209 .iter()
210 .map(|bands| {
211 bands
212 .berry_curvature
213 .iter()
214 .zip(&bands.energies)
215 .map(|(&omega, &energy)| {
216 omega * occupation.value_unchecked(energy, mu)
217 })
218 .sum::<f64>()
219 })
220 .sum();
221 sum / nk as f64 / determinant
222 })
223 .collect();
224 Array1::from_vec(values)
225 }
226 Integration::EnergyCut => {
227 let chemical_potentials = Array1::from_iter(params.mu.iter().copied());
228 let kvec = crate::kpoints::gen_kmesh(&k_mesh)?;
229 let mut all_pts: Vec<VertexKernel> = (0..kvec.nrows())
230 .into_par_iter()
231 .map(|ik| {
232 self.compute_velocity_kernel(
233 &kvec.row(ik).to_owned(),
234 &dir_a,
235 &dir_b,
236 None,
237 Gauge::Atom,
238 spin,
239 )
240 })
241 .collect();
242 global_band_track(&mut all_pts, ¶ms.kmesh);
243 let width = occupation.energy_width()?;
244 let sigma = match DIM {
245 2 => integrate_fermi_cut_2d(
246 &all_pts,
247 &k_mesh,
248 &chemical_potentials,
249 width,
250 params.eta,
251 ),
252 3 => integrate_fermi_cut_3d(
253 &all_pts,
254 &k_mesh,
255 &chemical_potentials,
256 width,
257 params.eta,
258 ),
259 _ => {
260 return Err(crate::TbError::InvalidDimension {
261 dim: DIM,
262 supported: vec![2, 3],
263 });
264 }
265 };
266 sigma / determinant
267 }
268 Integration::Simplex => unreachable!("rejected during validation"),
269 };
270
271 Ok(HallConductivityResult {
272 chemical_potentials: params.mu.clone(),
273 conductivity,
274 })
275 }
276}