Struct Rustb::Model

pub struct Model {
    pub dim_r: usize,
    pub norb: usize,
    pub nsta: usize,
    pub natom: usize,
    pub spin: bool,
    pub lat: Array2<f64>,
    pub orb: Array2<f64>,
    pub atom: Array2<f64>,
    pub atom_list: Vec<usize>,
    pub ham: Array3<Complex<f64>>,
    pub hamR: Array2<isize>,
    pub rmatrix: Array4<Complex<f64>>,
}

这个是 tight-binding 模型的基本单位

Fields

dim_r: usize

• The real space dimension of the model.

norb: usize

• The number of orbitals in the model.

nsta: usize

• The number of states in the model. If spin is enabled, nsta=norb$\times$2

natom: usize

• The number of atoms in the model. The atom and atom_list at the back are used to store the positions of the atoms, and the number of orbitals corresponding to each atom.

spin: bool

• Whether the model has spin enabled. If enabled, spin=true

lat: Array2<f64>

• The lattice vector of the model, a dim_r$\times$dim_r matrix, the axis0 direction stores a 1$\times$dim_r lattice vector.

orb: Array2<f64>

• The position of the orbitals in the model. We use fractional coordinates uniformly.

atom: Array2<f64>

• The position of the atoms in the model, also in fractional coordinates.

atom_list: Vec<usize>

• The number of orbitals in the atoms, in the same order as the atom positions.

ham: Array3<Complex<f64>>

• The Hamiltonian of the model, $\bra{m0}\hat H\ket{nR}$, a three-dimensional complex tensor of size n_R$\times$nsta$\times$ nsta, where the first nsta*nsta matrix corresponds to hopping within the unit cell, i.e. <m0|H|n0>, and the subsequent matrices correspond to hopping within hamR.

hamR: Array2<isize>

• The distance between the unit cell hoppings, i.e. R in $\bra{m0}\hat H\ket{nR}$.

rmatrix: Array4<Complex<f64>>

• The position matrix, i.e. $\bra{m0}\hat{\bm r}\ket{nR}$.

Implementations

impl Model

这个 impl 是给 tight-binding 模型提供基础的函数.

pub fn tb_model( dim_r: usize, lat: Array2<f64>, orb: Array2<f64>, spin: bool, atom: Option<Array2<f64>>, atom_list: Option<Vec<usize>> ) -> Model

This function is used to initialize a Model. The variables that need to be input are as follows:

• dim_r: the dimension of the model

• lat: the lattice constant

• orb: the orbital coordinates

• spin: whether to consider spin

• atom: the atomic coordinates, which can be None

• atom_list: the number of orbitals for each atom, which can be None.

Note that if any of the atomic variables are None, it is better to make them all None.

pub fn set_hop<T: Data<Elem = isize>, U: hop_use>( &mut self, tmp: U, ind_i: usize, ind_j: usize, R: &ArrayBase<T, Ix1>, pauli: isize )

This function is used to add hopping to the model. The “set” indicates that it can be used to override previous hopping.

• tmp: the parameters for hopping

• ind_i and ind_j: the orbital indices in the Hamiltonian, representing hopping from i to j

• R: the position of the target unit cell for hopping

• pauli: can take the values of 0, 1, 2, or 3, representing $\sigma_0$, $\sigma_x$, $\sigma_y$, $\sigma_z$.

In general, this function is used to set $\bra{i\bm 0}\hat H\ket{j\bm R}=$tmp.

pub fn add_hop<T: Data<Elem = isize>, U: hop_use>( &mut self, tmp: U, ind_i: usize, ind_j: usize, R: &ArrayBase<T, Ix1>, pauli: isize )

参数和 set_hop 一致, 但是 $\bra{i\bm 0}\hat H\ket{j\bm R}$+=tmp

pub fn add_element( &mut self, tmp: Complex<f64>, ind_i: usize, ind_j: usize, R: &Array1<isize> )

参数和 set_hop 一致, 但是 $\bra{i\bm 0}\hat H\ket{j\bm R}$+=tmp

pub fn set_onsite(&mut self, tmp: &Array1<f64>, pauli: isize)

直接对对角项进行设置

pub fn add_onsite(&mut self, tmp: &Array1<f64>, pauli: isize)

直接对对角项进行设置

pub fn set_onsite_one(&mut self, tmp: f64, ind: usize, pauli: isize)

对 $\bra{i\bm 0}\hat H\ket{i\bm 0}$ 进行设置

pub fn del_hop( &mut self, ind_i: usize, ind_j: usize, R: &Array1<isize>, pauli: isize )

删除 $\bra{i\bm 0}\hat H\ket{j\bm R}$

pub fn k_path( &self, path: &Array2<f64>, nk: usize ) -> (Array2<f64>, Array1<f64>, Array1<f64>)

根据高对称点来生成高对称路径, 画能带图

pub fn gen_ham<S: Data<Elem = f64>>( &self, kvec: &ArrayBase<S, Ix1> ) -> Array2<Complex<f64>>

这个是做傅里叶变换, 将实空间的哈密顿量变换到倒空间的哈密顿量

具体来说, 就是

对于wannier 基函数 $\ket{\bm R,i}$, 我们和布洛赫函数 $\ket{\bm k,i}$ 的变换形式为 $$\ket{\bm k,i}=\sum_{\bm R} e^{i\bm k\cdot(\bm R+\bm\tau_i)}\ket{\bm R}$$ 这样, 我们有 $$\begin{aligned}H_{mn,\bm k}&=\bra{m\bm k}\hat H\ket{n\bm k}=\sum_{\bm R^\prime}\sum_{\bm R} \bra{m\bm R^\prime}\hat H\ket{n\bm R}e^{-i(\bm R’-\bm R+\bm\tau_m-\bm \tau_n)\cdot\bm k}\\ &=\sum_{\bm R} \bra{m\bm 0}\hat H\ket{n\bm R}e^{i(\bm R-\bm\tau_m+\bm \tau_n)\cdot\bm k} \end{aligned}$$

pub fn gen_r<S: Data<Elem = f64>>( &self, kvec: &ArrayBase<S, Ix1> ) -> Array3<Complex<f64>>

和 gen_ham 类似, 将 $\hat{\bm r}$ 进行傅里叶变换

使用这个函数需要wannier90 中设置 warite_rmn=true, 在拟合的时候得到 seedname_r.dat 文件.

$$\bm r_{mn,\bm k}=\bra{m\bm k}\hat{\bm r}\ket{n\bm k}=\sum_{\bm R} \bra{m\bm 0}\hat{\bm r}\ket{n\bm R}e^{i(\bm R-\bm\tau_i+\bm \tau_j)\cdot\bm k}$$

pub fn gen_v<S: Data<Elem = f64>>( &self, kvec: &ArrayBase<S, Ix1> ) -> Array3<Complex<f64>>

这个函数是用来生成速度算符的, 即 $\bra{m\bm k}\p_\ap H_{\bm k}\ket{n\bm k},$ 这里的基函数是布洛赫波函数

这里速度算符的计算公式, 我们在程序中采用 tight-binding 模型, 即傅里叶变换的时候考虑原子位置.

这样我们就有

$$ \begin{aligned} \bra{m\bm k}\p_\ap H_{\bm k}\ket{n\bm k}&=\p_\ap\left(\bra{m\bm k} H\ket{n\bm k}\rt)-\p_\ap\left(\bra{m\bm k}\rt) H\ket{n\bm k}-\bra{m\bm k} H\p_\ap\ket{n\bm k}\\ &=\sum_{\bm R} i(\bm R-\bm\tau_m+\bm\tau_n)H_{mn}(\bm R) e^{i\bm k\cdot(\bm R-\bm\tau_m+\bm\tau_n)}-\lt[H_{\bm k},\mathcal A_{\bm k,\ap}\rt]_{mn} \end{aligned} $$

这里的 $\mathcal A_{\bm k}$ 的定义为 $$\mathcal A_{\bm k,\ap,mn}=-i\sum_{\bm R}r_{mn,\ap}(\bm R)e^{i\bm k\cdot(\bm R+\bm\tau_m-\bm\tau_{n})}+i\tau_{n\ap}\dt_{mn}$$ 其中 $\bm r_{mn}$ 可以由 wannier90 给出, 只需要设定 ``write_rmn=ture“ 在这里, 所有的 $\bm R$, $\bm r$, 以及 $\bm \tau$ 都是以实空间为坐标系.

pub fn solve_band_onek<S: Data<Elem = f64>>( &self, kvec: &ArrayBase<S, Ix1> ) -> Array1<f64>

求解单个k点的能带值

pub fn solve_band_range_onek<S: Data<Elem = f64>>( &self, kvec: &ArrayBase<S, Ix1>, range: (f64, f64), epsilon: f64 ) -> Array1<f64>

pub fn solve_band_all<S: Data<Elem = f64>>( &self, kvec: &ArrayBase<S, Ix2> ) -> Array2<f64>

求解多个k点的能带值

pub fn solve_band_all_parallel<S: Data<Elem = f64>>( &self, kvec: &ArrayBase<S, Ix2> ) -> Array2<f64>

并行求解多个k点的能带值

pub fn solve_onek<S: Data<Elem = f64>>( &self, kvec: &ArrayBase<S, Ix1> ) -> (Array1<f64>, Array2<Complex<f64>>)

pub fn solve_range_onek<S: Data<Elem = f64>>( &self, kvec: &ArrayBase<S, Ix1>, range: (f64, f64), epsilon: f64 ) -> (Array1<f64>, Array2<Complex<f64>>)

pub fn solve_all<S: Data<Elem = f64>>( &self, kvec: &ArrayBase<S, Ix2> ) -> (Array2<f64>, Array3<Complex<f64>>)

pub fn solve_all_parallel<S: Data<Elem = f64>>( &self, kvec: &ArrayBase<S, Ix2> ) -> (Array2<f64>, Array3<Complex<f64>>)

pub fn cut_piece(&self, num: usize, dir: usize) -> Model

这个函数是用来将model的某个方向进行截断的

num:截出多少个原胞

dir:方向

返回一个model, 其中 dir 和输入的model是一致的, 但是轨道数目和原子数目都会扩大num倍, 沿着dir方向没有胞间hopping. This function is used to truncate a certain direction of a model.

Parameters:

• num: number of unit cells to truncate.
• dir: the direction to be truncated.

Returns a new model with the same direction as the input model, but with the number of orbitals and atoms increased by a factor of “num”. There is no inter-cell hopping along the “dir” direction.

pub fn cut_dot( &self, num: usize, shape: usize, dir: Option<Vec<usize>> ) -> Model

这个是用来且角态或者切棱态的

pub fn remove_orb(&mut self, orb_list: &Vec<usize>)

pub fn remove_atom(&mut self, atom_list: &Vec<usize>)

pub fn unfold( &self, U: &Array2<f64>, path: &Array2<f64>, nk: usize, E_min: f64, E_max: f64, E_n: usize, eta: f64, precision: f64 ) -> Array2<f64>

能带反折叠算法, 用来计算能带反折叠后的能带. 可以用来计算合金以及一些超胞 算法参考

首先, 我们定义超胞布里渊区下的哈密顿量 $H_{\bm K}$ 以及其格林函数 $$G(\og,\bm K)=(\og+i\eta-H_{\bm K})^{-1}$$

这里 $H_{\bm k}$ 是超胞的哈密顿量. 其本征值和本征态为 $\ve_{N\bm K}$ 和 $\bra{\psi_{N\bm K}}$

故我们可以在本征态下将格林函数写为 $$G(\og,\bm K)=\sum_{N}\f{\dyad{\psi_{N\bm K}}}{\og+i\eta-\ve_{N\bm K}}$$

再利用普函数定理, 有 $A(\og,\bm K)=-\f{1}{\pi}\Im G(\og,\bm K)$, 对其求trace, 我们就能画超胞的能谱.

但是, 我们希望得到的是原胞的能谱, 所以我们需要得到原胞的基, 即 $\ket{n\bm k}$.

反折叠后的能谱为 $$A_{nn}(\og,\bm k)=\sum_{N\bm K}\lt\vert \braket{n\bm k}{\psi_{N\bm K}}\rt\vert^2 A_{NN}(\og,\bm K)$$

接下来我们计算 $\braket{n\bm k}{\psi_{N\bm K}}$

首先, 我们有$$ \lt\{ \begin{aligned} \ket{N\bm K}&=\f{1}{\sqrt{V}}\sum_{\bm R}e^{-i\bm K\cdot(\bm R+\bm\tau_N)}\ket{N\bm R}\\ \ket{n\bm k}&=\f{1}{\sqrt{v}}\sum_{\bm r}e^{-i\bm k\cdot(\bm r+\bm\tau_n)}\ket{n\bm r}\\ \end{aligned}\rt.$$

pub fn make_supercell(&self, U: &Array2<f64>) -> Model

This function is used to transform the model, where the new basis after transformation is given by $L’ = UL$.

pub fn shift_to_zero(&mut self)

这个是把超出单元格的原子移动回原本的位置的代码 这个可以用来重新构造原胞原子的位置. 我们看一下 SSH model 的结果

pub fn dos( &self, k_mesh: &Array1<usize>, E_min: f64, E_max: f64, E_n: usize, sigma: f64 ) -> (Array1<f64>, Array1<f64>)

我这里用的算法是高斯算法, 其算法过程如下 首先, 根据 k_mesh 算出所有的能量 $\ve_n$, 然后, 按照定义 $$\rho(\ve)=\sum_N\int\dd\bm k \delta(\ve_n-\ve)$$ 我们将 $\delta(\ve_n-\ve)$ 做了替换, 换成了 $\f{1}{\sqrt{2\pi}\sigma}e^{-\f{(\ve_n-\ve)^2}{2\sigma^2}}$ 然后, 计算方法是先算出所有的能量, 再将能量乘以高斯分布, 就能得到态密度. 态密度的光滑程度和k点密度以及高斯分布的展宽有关

pub fn show_band( &self, path: &Array2<f64>, label: &Vec<&str>, nk: usize, name: &str ) -> Result<()>

这个函数是用来快速画能带图的, 用python画图, 因为Rust画图不太方便.

pub fn from_hr(path: &str, file_name: &str, zero_energy: f64) -> Model

这个函数是从 wannier90 中读取 TB 文件的

这里 path 表示 文件的位置, 可以使用绝对路径, 即 “/” 开头的路径, 也可以使用相对路径, 即运行 cargo run 时候的文件夹作为起始路径.

file_name 就是 wannier90 中的 seedname, 文件可以读取 seedname.win, seedname_centres.xyz, seedname_hr.dat, 以及可选的 seedname_r.dat.

这里 seedname_centres.xyz 需要在 wannier90 中设置 write_xyz=true, 而 seedname_hr.dat 需要设置 write_hr=true.

impl Model

这个模块是用来提供电导率张量的, 包括自旋霍尔电导率和霍尔电导率, 以及非线性霍尔电导率.

pub fn berry_curvature_n_onek<S: Data<Elem = f64>>( &self, k_vec: &ArrayBase<S, Ix1>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, og: f64, spin: usize, eta: f64 ) -> (Array1<f64>, Array1<f64>)

给定一个k点, 返回 $\Omega_n(\bm k)$

pub fn berry_curvature_onek<S: Data<Elem = f64>>( &self, k_vec: &ArrayBase<S, Ix1>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, mu: f64, T: f64, og: f64, spin: usize, eta: f64 ) -> f64

给定一个 k 点, 指定 dir_1=$\alpha$, dir_2=$\beta$, T 代表温度, og= $\og$, mu=$\mu$ 为费米能级, spin=0,1,2,3 为$\sg_0,\sg_x,\sg_y,\sg_z$, 当体系不存在自旋的时候无论如何输入spin都默认 spin=0 eta=$\eta$ 是一个小量 这个函数返回的是 $$ \sum_n f_n\Omega_{n,\ap\bt}^\gm(\bm k)=\sum_n \f{1}{e^{(\ve_{n\bm k}-\mu)/T/k_B}+1} \sum_{m=\not n}\f{J_{\ap,nm}^\gm v_{\bt,mn}}{(\ve_{n\bm k}-\ve_{m\bm k})^2-(\og+i\eta)^2}$$ 其中 $J_\ap^\gm=\{s_\gm,v_\ap\}$

pub fn berry_curvature<S: Data<Elem = f64>>( &self, k_vec: &ArrayBase<S, Ix2>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, mu: f64, T: f64, og: f64, spin: usize, eta: f64 ) -> Array1<f64>

这个是用来并行计算大量k点的贝利曲率 这个可以用来画能带上的贝利曲率, 或者画一个贝利曲率的热图

pub fn Hall_conductivity( &self, k_mesh: &Array1<usize>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, mu: f64, T: f64, og: f64, spin: usize, eta: f64 ) -> f64

这个是用来计算霍尔电导的. 这里采用的是均匀撒点的方法, 利用 berry_curvature, 我们有 $$\sg_{\ap\bt}^\gm=\f{1}{N(2\pi)^r V}\sum_{\bm k} \Og_{\ap\bt}(\bm k),$$ 其中 $N$ 是 k 点数目,

pub fn Hall_conductivity_adapted( &self, k_mesh: &Array1<usize>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, mu: f64, T: f64, og: f64, spin: usize, eta: f64, re_err: f64, ab_err: f64 ) -> f64

这个是采用自适应积分算法来计算霍尔电导的, 一般来说, 我们建议 re_err 设置为 1, 而 ab_err 设置为 0.01

pub fn Hall_conductivity_mu( &self, k_mesh: &Array1<usize>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, mu: &Array1<f64>, T: f64, og: f64, spin: usize, eta: f64 ) -> Array1<f64>

用来计算多个 $\mu$ 值的, 这个函数是先求出 $\Omega_n$, 然后再分别用不同的费米能级来求和, 这样速度更快, 因为避免了重复求解 $\Omega_n$, 但是相对来说更耗内存, 而且不能做到自适应积分算法.

pub fn berry_curvature_dipole_n_onek( &self, k_vec: &Array1<f64>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, dir_3: &Array1<f64>, og: f64, spin: usize, eta: f64 ) -> (Array1<f64>, Array1<f64>)

这个是用来计算 $$\pdv{\ve_{n\bm k}}{k_\gm}\Og_{n,\ap\bt}$$

这里需要注意的一点是, 一般来说对于 $\p_\ap\ve_{\bm k}$, 需要用差分法来求解, 我这里提供了一个算法. $$ \ve_{\bm k}=U^\dag H_{\bm k} U\Rightarrow \pdv{\ve_{\bm k}}{\bm k}=U^\dag\pdv{H_{\bm k}}{\bm k}U+\pdv{U^\dag}{\bm k} H_{\bm k}U+U^\dag H_{\bm k}\pdv{U}{\bm k}$$ 因为 $U^\dag U=1\Rightarrow \p_{\bm k}U^\dag U=-U^\dag\p_{\bm k}U$, $\p_{\bm k}H_{\bm k}=v_{\bm k}$我们有 $$\pdv{\ve_{\bm k}}{\bm k}=v_{\bm k}+\lt[\ve_{\bm k},U^\dag\p_{\bm k}U\rt]$$ 而这里面唯一比较难求的项是 $D_{\bm k}=U^\dag\p_{\bm k}U$. 按照 vanderbilt 2008 年的论文中的公式, 用微扰论有 $$D_{mn,\bm k}=\left\{\begin{aligned}\f{v_{mn,\bm k}}{\ve_n-\ve_m} \quad &\text{if}\ m\ =\not n\\ 0 \quad \quad &\text{if}\ m\ = n\end{aligned}\right.$$ 我们观察到第二项对对角部分没有贡献, 所以我们可以直接设置为 $$\pdv{\ve_{\bm k}}{\bm k}=\text{diag}\lt(v_{\bm k}\rt)$$

pub fn berry_curvature_dipole_n( &self, k_vec: &Array2<f64>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, dir_3: &Array1<f64>, og: f64, spin: usize, eta: f64 ) -> (Array2<f64>, Array2<f64>)

This function performs parallel computation based on the onek function to obtain a series of Berry curvature dipoles at different k-points. 这个方法用的是对费米分布的修正, 因为高阶的dipole 修正导致的非线性霍尔电导为 $$\sg_{\ap\bt\gm}=\tau\int\dd\bm k\sum_n\p_\gm\ve_{n\bm k}\Og_{n,\ap\bt}\lt.\pdv{f_{\bm k}}{\ve}\rt\rvert_{E=\ve_{n\bm k}}.$$ 所以我们这里输出的是 $$\p_\gm\ve_{n\bm k}\Og_{n,\ap\bt}.$$

pub fn Nonlinear_Hall_conductivity_Extrinsic( &self, k_mesh: &Array1<usize>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, dir_3: &Array1<f64>, mu: &Array1<f64>, T: f64, og: f64, spin: usize, eta: f64 ) -> Array1<f64>

This function calculates the extrinsic nonlinear Hall conductivity by integrating the Berry curvature dipole over the entire Brillouin zone. The Berry curvature dipole is first computed at a series of k-points using parallel computation based on the onek function. 我们基于 berry_curvature_n_dipole 来并行得到所有 k 点的 $\p_\gm\ve_{n\bm k}\Og_{n,\ap\bt}$, 但是我们最后的公式为 $$\mathcal D_{\ap\bt\gm}=\int \dd\bm k \sum_n\lt(-\pdv{f_{n}}{\ve}\rt)\p_\gm\ve_{n\bm k}\Og_{n,\ap\bt}$$ 然而, $$-\pdv{f_{n}}{\ve}=\beta\f{e^{beta(\ve_n-\mu)}}{(e^{beta(\ve_n-\mu)}+1)^2}=\beta f_n(1-f_n)$$ 对于 T=0 的情况, 我们将采用四面体积分来替代, 这个需要很高的k点密度, 不建议使用 对于 T!=0 的情况, 我们会采用类似 Dos 的方法来计算

pub fn berry_connection_dipole_onek( &self, k_vec: &Array1<f64>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, dir_3: &Array1<f64>, spin: usize ) -> (Array1<f64>, Array1<f64>, Option<Array1<f64>>)

这个是根据 Nonlinear_Hall_conductivity_intrinsic 的注释, 当不存在自旋的时候提供 $$v_\ap G_{\bt\gm}-v_\bt G_{\ap\gm}$$ 其中 $$ G_{ij}=-2\text{Re}\sum_{m=\not n}\f{v_{i,nm}v_{j,mn}}{\lt(\ve_n-\ve_m\rt)^3} $$ 如果存在自旋, 即spin不等于0, 则还存在 $\p_{h_i} G_{jk}$ 项, 具体请看下面的非线性霍尔部分 我们这里暂时不考虑磁场, 只考虑电场

pub fn berry_connection_dipole( &self, k_vec: &Array2<f64>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, dir_3: &Array1<f64>, spin: usize ) -> (Array2<f64>, Array2<f64>, Option<Array2<f64>>)

这个是基于 onek 的, 进行关于 k 点并行求解

pub fn Nonlinear_Hall_conductivity_Intrinsic( &self, k_mesh: &Array1<usize>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, dir_3: &Array1<f64>, mu: &Array1<f64>, T: f64, spin: usize ) -> Array1<f64>

The Intrinsic Nonlinear Hall Conductivity arises from the correction of the Berry connection by the electric and magnetic fields [PRL 112, 166601 (2014)]. The formula employed is: $$\tilde\bm\Og_{\bm k}=\nb_{\bm k}\times\lt(\bm A_{\bm k}+\bm A_{\bm k}^\prime\rt)$$ and the $\bm A_{i,\bm k}^\prime=F_{ij}B_j+G_{ij}E_j$, where $$ \begin{aligned} F_{ij}&=\text{Im}\sum_{m=\not n}\f{v_{i,nm}\og_{j,mn}}{\lt(\ve_{n}-\ve_m\rt)^2}\\ G_{ij}&=2\text{Re}\sum_{m=\not n}\f{v_{i,nm}v_{j,mn}}{\lt(\ve_n-\ve_m\rt)^3}\\ \og_{\ap,mn}&=-i\ep_{\ap\bt\gm}\sum_{l=\not n}\f{\lt(v_{\bt,ml}+\p_\bt \ve_{\bm k}\dt_{ml}\rt)v_{\gm,ln}}{\ve_l-\ve_n} \end{aligned} $$ 最后我们有 $$ \bm j^\prime=\bm E\times\int\f{\dd\bm k}{(2\pi)^3}\lt[\p_{\bm k}\ve_{\bm k}\times\bm A^\prime+\bm\Og\lt(\bm B\cdot\bm m\rt)\rt]\pdv{f_{\bm k}}{\ve} $$ 对其对电场和磁场进行偏导, 有 $$ \begin{aligned} \f{\p^2 j_{\ap}^\prime}{\p E_\bt\p E_\gm}&=\int\f{\dd\bm k}{(2\pi)^3}\lt(\p_\ap\ve_{\bm k} G_{\bt\gm}-\p_\bt\ve_{\bm k} G_{\ap\gm}\rt)\pdv{f_{\bm k}}{\ve}\\ \f{\p^2 j_{\ap}^\prime}{\p E_\bt\p B_\gm}&=\int\f{\dd\bm k}{(2\pi)^3}\lt(\p_\ap\ve_{\bm k} F_{\bt\gm}-\p_\bt\ve_{\bm k} F_{\ap\gm}+\ep_{\ap\bt\ell}\Og_{\ell} m_\gm\rt)\pdv{f_{\bm k}}{\ve} \end{aligned} $$ 由于存在 $\pdv{f_{\bm k}}{\ve}$, 不建议将温度 T=0

可以考虑当 T=0 时候, 利用高斯公式, 将费米面内的部分进行积分, 得到精确解. 但是我现在还没办法很好的求解费米面, 所以暂时不考虑这个算法.而且对于二维体系, 公式还不一样, 还得分步讨论, 后面有时间再考虑这个程序.

对于自旋霍尔效应, 按照文章 [PRL 112, 166601 (2014)], 非线性自旋霍尔电导为 $$\sg_{\ap\bt\gm}^i=-\int\dd\bm k \lt[\f{1}{2}f_{\bm k}\pdv{G_{\bt\gm}}{h_\ap}+\pdv{f_{\bm k}}{\ve}\lt(\p_{\ap}s_{\bm k}^i G_{\bt\gm}-\p_\bt\ve_{\bm k}G_{\ap\gm}^h\rt)\rt]$$ 其中 $$\f{\p G_{\bt\gm,n}}{\p h_\ap}=2\text{Re}\sum_{n^\pr =\not n}\f{3\lt(s^i_{\ap,n}-s^i_{\ap,n_1}\rt)v_{\bt,nn_1} v_{\gm,n^\pr n}}{\lt(\ve_n-\ve_{n^\pr}\rt)^4}-2\text{Re}\sum_{n_1=\not n}\sum_{n_2=\not n}\lt[\f{s^i_{\ap,nn_2} v_{\bt,n_2n_1} v_{\gm,n_1 n}}{\lt(\ve_n-\ve_{n_1}\rt)^3(\ve_n-\ve_{n_2})}+(\bt \leftrightarrow \gm)\rt]-2\text{Re}\sum_{n_1=\not n}\sum_{n_2=\not n_1}\lt[\f{s^i_{\ap,n_1n_2} v_{\bt,n_2n} v_{\gm,n n_1}}{\lt(\ve_n-\ve_{n_1}\rt)^3(\ve_{n_1}-\ve_{n_2})}+(\bt \leftrightarrow \gm)\rt]$$ 以及 $$ \lt\{\begin{aligned} G_{\ap\bt}&=2\text{Re}\sum_{m=\not n}\f{v_{\ap,nm}v_{\bt,mn}}{\lt(\ve_n-\ve_m\rt)^3}\\ G_{\ap\bt}^h&=2\text{Re}\sum_{m=\not n}\f{s^i_{\ap,nm}v_{\bt,mn}}{\lt(\ve_n-\ve_m\rt)^3}\\ \end{aligned}\rt. $$

这里 $s^i_{\ap,mn}$ 的具体形式, 原文中没有明确给出, 但是我根据霍尔效应的类比, 我猜是 $\{\hat s^i,v_\ap\}$

impl Model

这个模块是用wilson loop 的方法来计算各种几何量.

pub fn berry_loop<S>( &self, kvec: &ArrayBase<S, Ix2>, occ: &Vec<usize> ) -> Array1<f64>

where S: Data<Elem = f64>,

这个函数是计算某一个闭合路径上的 berry phase, 用的是wilson loop 方法.l

其算法如下: 我们首先将末端的k点的波函数相位统一. 如果闭合回路沿着布里渊区两端, 那么因为布洛赫函数在布里渊区存在一个相位差 $e^{-2\pi i \bm \tau_i}$, 我们有 $\ket{\psi_{n,\bm k_\text{end}}}=e^{-2\pi i \bm \tau_i}\ket{\psi_{n,\bm k_\text{first}}}$

我们定义交叠矩阵 $F_{mn,\bm k}=\braket{\psi_{m,\bm k}}{\psi_{n,\bm k+\dd\bm k}}$

接下来我们将其正交化, 用 SVD 分解, 有 $U,S,V=\text{svd}(F_{\bm k})$, $F_{\bm k}=UV$

接下来我们 将其 连乘, 有 $$W=\prod_{i} F_{\bm k_i} F_{\bm k_{i+1}}$$,

最后, 我们求本征值并取其幅角, 有 $e^{i\Theta}=\text{eigh}(W)$, 就能够得到这个loop的wanniercentre

pub fn berry_flux( &self, occ: &Vec<usize>, k_start: &Array1<f64>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, nk1: usize, nk2: usize ) -> Array3<f64>

这个函数是用 wilson loop 方法来计算berry curvature 的. 根据给定的平面, 其返回一个 Array2, 这个算法的优点是精度高, 计算量小, 能快速收敛, 但是只能用于绝缘体. 前两个指标表示横和纵, 数值表示大小

pub fn berry_phase(&self, occ: &Vec<usize>, kvec: &Array3<f64>) -> Array2<f64>

这里是计算wcc的, 沿着第一个方向走, 沿着第二个方向积分

pub fn wannier_centre( &self, occ: &Vec<usize>, k_start: &Array1<f64>, dir_1: &Array1<f64>, dir_2: &Array1<f64>, nk1: usize, nk2: usize ) -> Array2<f64>

这里是计算wcc的, 沿着第一个方向走, 沿着第二个方向积分

Trait Implementations

impl Clone for Model

fn clone(&self) -> Model

Returns a copy of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Debug for Model

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.