diff --git a/PyPI/README.md b/PyPI/README.md index c6f4673..19e90a2 100644 --- a/PyPI/README.md +++ b/PyPI/README.md @@ -1 +1 @@ -Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about (Ji-Huan Guan, 关济寰). The primary location of this package is on website https://py.guanjihuan.com. \ No newline at end of file +Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about (Ji-Huan Guan, 关济寰). The primary location of this package is on website https://py.guanjihuan.com. The GitHub location of this package is on https://github.com/guanjihuan/py.guanjihuan.com. \ No newline at end of file diff --git a/PyPI/setup.cfg b/PyPI/setup.cfg index a4a91f3..06cf1da 100644 --- a/PyPI/setup.cfg +++ b/PyPI/setup.cfg @@ -1,7 +1,7 @@ [metadata] # replace with your username: name = guan -version = 0.1.12 +version = 0.1.13 author = guanjihuan author_email = guanjihuan@163.com description = An open source python package diff --git a/PyPI/src/guan.egg-info/PKG-INFO b/PyPI/src/guan.egg-info/PKG-INFO index 138a58c..1ba1fb9 100644 --- a/PyPI/src/guan.egg-info/PKG-INFO +++ b/PyPI/src/guan.egg-info/PKG-INFO @@ -1,6 +1,6 @@ Metadata-Version: 2.1 Name: guan -Version: 0.1.12 +Version: 0.1.13 Summary: An open source python package Home-page: https://py.guanjihuan.com Author: guanjihuan @@ -13,4 +13,4 @@ Requires-Python: >=3.6 Description-Content-Type: text/markdown License-File: LICENSE -Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about (Ji-Huan Guan, 关济寰). The primary location of this package is on website https://py.guanjihuan.com. +Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about (Ji-Huan Guan, 关济寰). The primary location of this package is on website https://py.guanjihuan.com. The GitHub location of this package is on https://github.com/guanjihuan/py.guanjihuan.com. diff --git a/PyPI/src/guan.egg-info/SOURCES.txt b/PyPI/src/guan.egg-info/SOURCES.txt index d26d474..6a00c55 100644 --- a/PyPI/src/guan.egg-info/SOURCES.txt +++ b/PyPI/src/guan.egg-info/SOURCES.txt @@ -2,7 +2,20 @@ LICENSE README.md pyproject.toml setup.cfg +src/guan/Fourier_transform.py +src/guan/Green_functions.py +src/guan/Hamiltonian_of_finite_size_systems.py +src/guan/Hamiltonian_of_models_in_reciprocal_space.py src/guan/__init__.py +src/guan/band_structures_and_wave_functions.py +src/guan/basic_functions.py +src/guan/data_processing.py +src/guan/density_of_states.py +src/guan/file_processing.py +src/guan/plot_figures.py +src/guan/quantum_transport.py +src/guan/read_and_write.py +src/guan/topological_invariant.py src/guan.egg-info/PKG-INFO src/guan.egg-info/SOURCES.txt src/guan.egg-info/dependency_links.txt diff --git a/PyPI/src/guan/Fourier_transform.py b/PyPI/src/guan/Fourier_transform.py new file mode 100644 index 0000000..50912f4 --- /dev/null +++ b/PyPI/src/guan/Fourier_transform.py @@ -0,0 +1,135 @@ +# Module: Fourier_transform + +# 通过元胞和跃迁项得到一维的哈密顿量（需要输入k值） +def one_dimensional_fourier_transform(k, unit_cell, hopping): + import numpy as np + import cmath + unit_cell = np.array(unit_cell) + hopping = np.array(hopping) + hamiltonian = unit_cell+hopping*cmath.exp(1j*k)+hopping.transpose().conj()*cmath.exp(-1j*k) + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 通过元胞和跃迁项得到二维方格子的哈密顿量（需要输入k值） +def two_dimensional_fourier_transform_for_square_lattice(k1, k2, unit_cell, hopping_1, hopping_2): + import numpy as np + import cmath + unit_cell = np.array(unit_cell) + hopping_1 = np.array(hopping_1) + hopping_2 = np.array(hopping_2) + hamiltonian = unit_cell+hopping_1*cmath.exp(1j*k1)+hopping_1.transpose().conj()*cmath.exp(-1j*k1)+hopping_2*cmath.exp(1j*k2)+hopping_2.transpose().conj()*cmath.exp(-1j*k2) + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 通过元胞和跃迁项得到三维立方格子的哈密顿量（需要输入k值） +def three_dimensional_fourier_transform_for_cubic_lattice(k1, k2, k3, unit_cell, hopping_1, hopping_2, hopping_3): + import numpy as np + import cmath + unit_cell = np.array(unit_cell) + hopping_1 = np.array(hopping_1) + hopping_2 = np.array(hopping_2) + hopping_3 = np.array(hopping_3) + hamiltonian = unit_cell+hopping_1*cmath.exp(1j*k1)+hopping_1.transpose().conj()*cmath.exp(-1j*k1)+hopping_2*cmath.exp(1j*k2)+hopping_2.transpose().conj()*cmath.exp(-1j*k2)+hopping_3*cmath.exp(1j*k3)+hopping_3.transpose().conj()*cmath.exp(-1j*k3) + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 通过元胞和跃迁项得到一维的哈密顿量（返回的哈密顿量为携带k的函数） +def one_dimensional_fourier_transform_with_k(unit_cell, hopping): + import functools + import guan + hamiltonian_function = functools.partial(guan.one_dimensional_fourier_transform, unit_cell=unit_cell, hopping=hopping) + guan.statistics_of_guan_package() + return hamiltonian_function + +# 通过元胞和跃迁项得到二维方格子的哈密顿量（返回的哈密顿量为携带k的函数） +def two_dimensional_fourier_transform_for_square_lattice_with_k1_k2(unit_cell, hopping_1, hopping_2): + import functools + import guan + hamiltonian_function = functools.partial(guan.two_dimensional_fourier_transform_for_square_lattice, unit_cell=unit_cell, hopping_1=hopping_1, hopping_2=hopping_2) + guan.statistics_of_guan_package() + return hamiltonian_function + +# 通过元胞和跃迁项得到三维立方格子的哈密顿量（返回的哈密顿量为携带k的函数） +def three_dimensional_fourier_transform_for_cubic_lattice_with_k1_k2_k3(unit_cell, hopping_1, hopping_2, hopping_3): + import functools + import guan + hamiltonian_function = functools.partial(guan.three_dimensional_fourier_transform_for_cubic_lattice, unit_cell=unit_cell, hopping_1=hopping_1, hopping_2=hopping_2, hopping_3=hopping_3) + guan.statistics_of_guan_package() + return hamiltonian_function + +# 由实空间格矢得到倒空间格矢（一维） +def calculate_one_dimensional_reciprocal_lattice_vector(a1): + import numpy as np + b1 = 2*np.pi/a1 + import guan + guan.statistics_of_guan_package() + return b1 + +# 由实空间格矢得到倒空间格矢（二维） +def calculate_two_dimensional_reciprocal_lattice_vectors(a1, a2): + import numpy as np + a1 = np.array(a1) + a2 = np.array(a2) + a1 = np.append(a1, 0) + a2 = np.append(a2, 0) + a3 = np.array([0, 0, 1]) + b1 = 2*np.pi*np.cross(a2, a3)/np.dot(a1, np.cross(a2, a3)) + b2 = 2*np.pi*np.cross(a3, a1)/np.dot(a1, np.cross(a2, a3)) + b1 = np.delete(b1, 2) + b2 = np.delete(b2, 2) + import guan + guan.statistics_of_guan_package() + return b1, b2 + +# 由实空间格矢得到倒空间格矢（三维） +def calculate_three_dimensional_reciprocal_lattice_vectors(a1, a2, a3): + import numpy as np + a1 = np.array(a1) + a2 = np.array(a2) + a3 = np.array(a3) + b1 = 2*np.pi*np.cross(a2, a3)/np.dot(a1, np.cross(a2, a3)) + b2 = 2*np.pi*np.cross(a3, a1)/np.dot(a1, np.cross(a2, a3)) + b3 = 2*np.pi*np.cross(a1, a2)/np.dot(a1, np.cross(a2, a3)) + import guan + guan.statistics_of_guan_package() + return b1, b2, b3 + +# 由实空间格矢得到倒空间格矢（一维），这里为符号运算 +def calculate_one_dimensional_reciprocal_lattice_vector_with_sympy(a1): + import sympy + b1 = 2*sympy.pi/a1 + import guan + guan.statistics_of_guan_package() + return b1 + +# 由实空间格矢得到倒空间格矢（二维），这里为符号运算 +def calculate_two_dimensional_reciprocal_lattice_vectors_with_sympy(a1, a2): + import sympy + a1 = sympy.Matrix(1, 3, [a1[0], a1[1], 0]) + a2 = sympy.Matrix(1, 3, [a2[0], a2[1], 0]) + a3 = sympy.Matrix(1, 3, [0, 0, 1]) + cross_a2_a3 = a2.cross(a3) + cross_a3_a1 = a3.cross(a1) + b1 = 2*sympy.pi*cross_a2_a3/a1.dot(cross_a2_a3) + b2 = 2*sympy.pi*cross_a3_a1/a1.dot(cross_a2_a3) + b1 = sympy.Matrix(1, 2, [b1[0], b1[1]]) + b2 = sympy.Matrix(1, 2, [b2[0], b2[1]]) + import guan + guan.statistics_of_guan_package() + return b1, b2 + +# 由实空间格矢得到倒空间格矢（三维），这里为符号运算 +def calculate_three_dimensional_reciprocal_lattice_vectors_with_sympy(a1, a2, a3): + import sympy + cross_a2_a3 = a2.cross(a3) + cross_a3_a1 = a3.cross(a1) + cross_a1_a2 = a1.cross(a2) + b1 = 2*sympy.pi*cross_a2_a3/a1.dot(cross_a2_a3) + b2 = 2*sympy.pi*cross_a3_a1/a1.dot(cross_a2_a3) + b3 = 2*sympy.pi*cross_a1_a2/a1.dot(cross_a2_a3) + import guan + guan.statistics_of_guan_package() + return b1, b2, b3 diff --git a/PyPI/src/guan/Green_functions.py b/PyPI/src/guan/Green_functions.py new file mode 100644 index 0000000..7f4713f --- /dev/null +++ b/PyPI/src/guan/Green_functions.py @@ -0,0 +1,153 @@ +# Module: Green_functions + +# 输入哈密顿量，得到格林函数 +def green_function(fermi_energy, hamiltonian, broadening, self_energy=0): + import numpy as np + if np.array(hamiltonian).shape==(): + dim = 1 + else: + dim = np.array(hamiltonian).shape[0] + green = np.linalg.inv((fermi_energy+broadening*1j)*np.eye(dim)-hamiltonian-self_energy) + import guan + guan.statistics_of_guan_package() + return green + +# 在Dyson方程中的一个中间格林函数G_{nn}^{n} +def green_function_nn_n(fermi_energy, h00, h01, green_nn_n_minus, broadening, self_energy=0): + import numpy as np + h01 = np.array(h01) + if np.array(h00).shape==(): + dim = 1 + else: + dim = np.array(h00).shape[0] + green_nn_n = np.linalg.inv((fermi_energy+broadening*1j)*np.identity(dim)-h00-np.dot(np.dot(h01.transpose().conj(), green_nn_n_minus), h01)-self_energy) + import guan + guan.statistics_of_guan_package() + return green_nn_n + +# 在Dyson方程中的一个中间格林函数G_{in}^{n} +def green_function_in_n(green_in_n_minus, h01, green_nn_n): + import numpy as np + green_in_n = np.dot(np.dot(green_in_n_minus, h01), green_nn_n) + import guan + guan.statistics_of_guan_package() + return green_in_n + +# 在Dyson方程中的一个中间格林函数G_{ni}^{n} +def green_function_ni_n(green_nn_n, h01, green_ni_n_minus): + import numpy as np + h01 = np.array(h01) + green_ni_n = np.dot(np.dot(green_nn_n, h01.transpose().conj()), green_ni_n_minus) + import guan + guan.statistics_of_guan_package() + return green_ni_n + +# 在Dyson方程中的一个中间格林函数G_{ii}^{n} +def green_function_ii_n(green_ii_n_minus, green_in_n_minus, h01, green_nn_n, green_ni_n_minus): + import numpy as np + green_ii_n = green_ii_n_minus+np.dot(np.dot(np.dot(np.dot(green_in_n_minus, h01), green_nn_n), h01.transpose().conj()),green_ni_n_minus) + import guan + guan.statistics_of_guan_package() + return green_ii_n + +# 计算转移矩阵（该矩阵可以用来计算表面格林函数） +def transfer_matrix(fermi_energy, h00, h01): + import numpy as np + h01 = np.array(h01) + if np.array(h00).shape==(): + dim = 1 + else: + dim = np.array(h00).shape[0] + transfer = np.zeros((2*dim, 2*dim), dtype=complex) + transfer[0:dim, 0:dim] = np.dot(np.linalg.inv(h01), fermi_energy*np.identity(dim)-h00) + transfer[0:dim, dim:2*dim] = np.dot(-1*np.linalg.inv(h01), h01.transpose().conj()) + transfer[dim:2*dim, 0:dim] = np.identity(dim) + transfer[dim:2*dim, dim:2*dim] = 0 + import guan + guan.statistics_of_guan_package() + return transfer + +# 计算电极的表面格林函数 +def surface_green_function_of_lead(fermi_energy, h00, h01): + import numpy as np + h01 = np.array(h01) + if np.array(h00).shape==(): + dim = 1 + else: + dim = np.array(h00).shape[0] + fermi_energy = fermi_energy+1e-9*1j + transfer = transfer_matrix(fermi_energy, h00, h01) + eigenvalue, eigenvector = np.linalg.eig(transfer) + ind = np.argsort(np.abs(eigenvalue)) + temp = np.zeros((2*dim, 2*dim), dtype=complex) + i0 = 0 + for ind0 in ind: + temp[:, i0] = eigenvector[:, ind0] + i0 += 1 + s1 = temp[dim:2*dim, 0:dim] + s2 = temp[0:dim, 0:dim] + s3 = temp[dim:2*dim, dim:2*dim] + s4 = temp[0:dim, dim:2*dim] + right_lead_surface = np.linalg.inv(fermi_energy*np.identity(dim)-h00-np.dot(np.dot(h01, s2), np.linalg.inv(s1))) + left_lead_surface = np.linalg.inv(fermi_energy*np.identity(dim)-h00-np.dot(np.dot(h01.transpose().conj(), s3), np.linalg.inv(s4))) + import guan + guan.statistics_of_guan_package() + return right_lead_surface, left_lead_surface + +# 计算电极的自能（基于Dyson方程的小矩阵形式） +def self_energy_of_lead(fermi_energy, h00, h01): + import numpy as np + import guan + h01 = np.array(h01) + right_lead_surface, left_lead_surface = guan.surface_green_function_of_lead(fermi_energy, h00, h01) + right_self_energy = np.dot(np.dot(h01, right_lead_surface), h01.transpose().conj()) + left_self_energy = np.dot(np.dot(h01.transpose().conj(), left_lead_surface), h01) + gamma_right = (right_self_energy - right_self_energy.transpose().conj())*1j + gamma_left = (left_self_energy - left_self_energy.transpose().conj())*1j + guan.statistics_of_guan_package() + return right_self_energy, left_self_energy, gamma_right, gamma_left + +# 计算电极的自能（基于中心区整体的大矩阵形式） +def self_energy_of_lead_with_h_LC_and_h_CR(fermi_energy, h00, h01, h_LC, h_CR): + import numpy as np + import guan + h_LC = np.array(h_LC) + h_CR = np.array(h_CR) + right_lead_surface, left_lead_surface = guan.surface_green_function_of_lead(fermi_energy, h00, h01) + right_self_energy = np.dot(np.dot(h_CR, right_lead_surface), h_CR.transpose().conj()) + left_self_energy = np.dot(np.dot(h_LC.transpose().conj(), left_lead_surface), h_LC) + gamma_right = (right_self_energy - right_self_energy.transpose().conj())*1j + gamma_left = (left_self_energy - left_self_energy.transpose().conj())*1j + guan.statistics_of_guan_package() + return right_self_energy, left_self_energy, gamma_right, gamma_left + +# 计算电极的自能（基于中心区整体的大矩阵形式，可适用于多端电导的计算） +def self_energy_of_lead_with_h_lead_to_center(fermi_energy, h00, h01, h_lead_to_center): + import numpy as np + import guan + h_lead_to_center = np.array(h_lead_to_center) + right_lead_surface, left_lead_surface = guan.surface_green_function_of_lead(fermi_energy, h00, h01) + self_energy = np.dot(np.dot(h_lead_to_center.transpose().conj(), right_lead_surface), h_lead_to_center) + gamma = (self_energy - self_energy.transpose().conj())*1j + guan.statistics_of_guan_package() + return self_energy, gamma + +# 计算考虑电极自能后的中心区的格林函数 +def green_function_with_leads(fermi_energy, h00, h01, h_LC, h_CR, center_hamiltonian): + import numpy as np + import guan + dim = np.array(center_hamiltonian).shape[0] + right_self_energy, left_self_energy, gamma_right, gamma_left = guan.self_energy_of_lead_with_h_LC_and_h_CR(fermi_energy, h00, h01, h_LC, h_CR) + green = np.linalg.inv(fermi_energy*np.identity(dim)-center_hamiltonian-left_self_energy-right_self_energy) + guan.statistics_of_guan_package() + return green, gamma_right, gamma_left + +# 计算用于计算局域电流的格林函数G_n +def electron_correlation_function_green_n_for_local_current(fermi_energy, h00, h01, h_LC, h_CR, center_hamiltonian): + import numpy as np + import guan + right_self_energy, left_self_energy, gamma_right, gamma_left = guan.self_energy_of_lead_with_h_LC_and_h_CR(fermi_energy, h00, h01, h_LC, h_CR) + green = guan.green_function(fermi_energy, center_hamiltonian, broadening=0, self_energy=left_self_energy+right_self_energy) + G_n = np.imag(np.dot(np.dot(green, gamma_left), green.transpose().conj())) + guan.statistics_of_guan_package() + return G_n \ No newline at end of file diff --git a/PyPI/src/guan/Hamiltonian_of_finite_size_systems.py b/PyPI/src/guan/Hamiltonian_of_finite_size_systems.py new file mode 100644 index 0000000..a8ffaa8 --- /dev/null +++ b/PyPI/src/guan/Hamiltonian_of_finite_size_systems.py @@ -0,0 +1,260 @@ +# Module: Hamiltonian_of_finite_size_systems + +# 构建一维的有限尺寸体系哈密顿量（可设置是否为周期边界条件） +def hamiltonian_of_finite_size_system_along_one_direction(N, on_site=0, hopping=1, period=0): + import numpy as np + on_site = np.array(on_site) + hopping = np.array(hopping) + if on_site.shape==(): + dim = 1 + else: + dim = on_site.shape[0] + hamiltonian = np.zeros((N*dim, N*dim), dtype=complex) + for i0 in range(N): + hamiltonian[i0*dim+0:i0*dim+dim, i0*dim+0:i0*dim+dim] = on_site + for i0 in range(N-1): + hamiltonian[i0*dim+0:i0*dim+dim, (i0+1)*dim+0:(i0+1)*dim+dim] = hopping + hamiltonian[(i0+1)*dim+0:(i0+1)*dim+dim, i0*dim+0:i0*dim+dim] = hopping.transpose().conj() + if period == 1: + hamiltonian[(N-1)*dim+0:(N-1)*dim+dim, 0:dim] = hopping + hamiltonian[0:dim, (N-1)*dim+0:(N-1)*dim+dim] = hopping.transpose().conj() + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 构建二维的方格子有限尺寸体系哈密顿量（可设置是否为周期边界条件） +def hamiltonian_of_finite_size_system_along_two_directions_for_square_lattice(N1, N2, on_site=0, hopping_1=1, hopping_2=1, period_1=0, period_2=0): + import numpy as np + on_site = np.array(on_site) + hopping_1 = np.array(hopping_1) + hopping_2 = np.array(hopping_2) + if on_site.shape==(): + dim = 1 + else: + dim = on_site.shape[0] + hamiltonian = np.zeros((N1*N2*dim, N1*N2*dim), dtype=complex) + for i1 in range(N1): + for i2 in range(N2): + hamiltonian[i1*N2*dim+i2*dim+0:i1*N2*dim+i2*dim+dim, i1*N2*dim+i2*dim+0:i1*N2*dim+i2*dim+dim] = on_site + for i1 in range(N1-1): + for i2 in range(N2): + hamiltonian[i1*N2*dim+i2*dim+0:i1*N2*dim+i2*dim+dim, (i1+1)*N2*dim+i2*dim+0:(i1+1)*N2*dim+i2*dim+dim] = hopping_1 + hamiltonian[(i1+1)*N2*dim+i2*dim+0:(i1+1)*N2*dim+i2*dim+dim, i1*N2*dim+i2*dim+0:i1*N2*dim+i2*dim+dim] = hopping_1.transpose().conj() + for i1 in range(N1): + for i2 in range(N2-1): + hamiltonian[i1*N2*dim+i2*dim+0:i1*N2*dim+i2*dim+dim, i1*N2*dim+(i2+1)*dim+0:i1*N2*dim+(i2+1)*dim+dim] = hopping_2 + hamiltonian[i1*N2*dim+(i2+1)*dim+0:i1*N2*dim+(i2+1)*dim+dim, i1*N2*dim+i2*dim+0:i1*N2*dim+i2*dim+dim] = hopping_2.transpose().conj() + if period_1 == 1: + for i2 in range(N2): + hamiltonian[(N1-1)*N2*dim+i2*dim+0:(N1-1)*N2*dim+i2*dim+dim, i2*dim+0:i2*dim+dim] = hopping_1 + hamiltonian[i2*dim+0:i2*dim+dim, (N1-1)*N2*dim+i2*dim+0:(N1-1)*N2*dim+i2*dim+dim] = hopping_1.transpose().conj() + if period_2 == 1: + for i1 in range(N1): + hamiltonian[i1*N2*dim+(N2-1)*dim+0:i1*N2*dim+(N2-1)*dim+dim, i1*N2*dim+0:i1*N2*dim+dim] = hopping_2 + hamiltonian[i1*N2*dim+0:i1*N2*dim+dim, i1*N2*dim+(N2-1)*dim+0:i1*N2*dim+(N2-1)*dim+dim] = hopping_2.transpose().conj() + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 构建三维的立方格子有限尺寸体系哈密顿量（可设置是否为周期边界条件） +def hamiltonian_of_finite_size_system_along_three_directions_for_cubic_lattice(N1, N2, N3, on_site=0, hopping_1=1, hopping_2=1, hopping_3=1, period_1=0, period_2=0, period_3=0): + import numpy as np + on_site = np.array(on_site) + hopping_1 = np.array(hopping_1) + hopping_2 = np.array(hopping_2) + hopping_3 = np.array(hopping_3) + if on_site.shape==(): + dim = 1 + else: + dim = on_site.shape[0] + hamiltonian = np.zeros((N1*N2*N3*dim, N1*N2*N3*dim), dtype=complex) + for i1 in range(N1): + for i2 in range(N2): + for i3 in range(N3): + hamiltonian[i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim, i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim] = on_site + for i1 in range(N1-1): + for i2 in range(N2): + for i3 in range(N3): + hamiltonian[i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim, (i1+1)*N2*N3*dim+i2*N3*dim+i3*dim+0:(i1+1)*N2*N3*dim+i2*N3*dim+i3*dim+dim] = hopping_1 + hamiltonian[(i1+1)*N2*N3*dim+i2*N3*dim+i3*dim+0:(i1+1)*N2*N3*dim+i2*N3*dim+i3*dim+dim, i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim] = hopping_1.transpose().conj() + for i1 in range(N1): + for i2 in range(N2-1): + for i3 in range(N3): + hamiltonian[i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim, i1*N2*N3*dim+(i2+1)*N3*dim+i3*dim+0:i1*N2*N3*dim+(i2+1)*N3*dim+i3*dim+dim] = hopping_2 + hamiltonian[i1*N2*N3*dim+(i2+1)*N3*dim+i3*dim+0:i1*N2*N3*dim+(i2+1)*N3*dim+i3*dim+dim, i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim] = hopping_2.transpose().conj() + for i1 in range(N1): + for i2 in range(N2): + for i3 in range(N3-1): + hamiltonian[i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim, i1*N2*N3*dim+i2*N3*dim+(i3+1)*dim+0:i1*N2*N3*dim+i2*N3*dim+(i3+1)*dim+dim] = hopping_3 + hamiltonian[i1*N2*N3*dim+i2*N3*dim+(i3+1)*dim+0:i1*N2*N3*dim+i2*N3*dim+(i3+1)*dim+dim, i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim] = hopping_3.transpose().conj() + if period_1 == 1: + for i2 in range(N2): + for i3 in range(N3): + hamiltonian[(N1-1)*N2*N3*dim+i2*N3*dim+i3*dim+0:(N1-1)*N2*N3*dim+i2*N3*dim+i3*dim+dim, i2*N3*dim+i3*dim+0:i2*N3*dim+i3*dim+dim] = hopping_1 + hamiltonian[i2*N3*dim+i3*dim+0:i2*N3*dim+i3*dim+dim, (N1-1)*N2*N3*dim+i2*N3*dim+i3*dim+0:(N1-1)*N2*N3*dim+i2*N3*dim+i3*dim+dim] = hopping_1.transpose().conj() + if period_2 == 1: + for i1 in range(N1): + for i3 in range(N3): + hamiltonian[i1*N2*N3*dim+(N2-1)*N3*dim+i3*dim+0:i1*N2*N3*dim+(N2-1)*N3*dim+i3*dim+dim, i1*N2*N3*dim+i3*dim+0:i1*N2*N3*dim+i3*dim+dim] = hopping_2 + hamiltonian[i1*N2*N3*dim+i3*dim+0:i1*N2*N3*dim+i3*dim+dim, i1*N2*N3*dim+(N2-1)*N3*dim+i3*dim+0:i1*N2*N3*dim+(N2-1)*N3*dim+i3*dim+dim] = hopping_2.transpose().conj() + if period_3 == 1: + for i1 in range(N1): + for i2 in range(N2): + hamiltonian[i1*N2*N3*dim+i2*N3*dim+(N3-1)*dim+0:i1*N2*N3*dim+i2*N3*dim+(N3-1)*dim+dim, i1*N2*N3*dim+i2*N3*dim+0:i1*N2*N3*dim+i2*N3*dim+dim] = hopping_3 + hamiltonian[i1*N2*N3*dim+i2*N3*dim+0:i1*N2*N3*dim+i2*N3*dim+dim, i1*N2*N3*dim+i2*N3*dim+(N3-1)*dim+0:i1*N2*N3*dim+i2*N3*dim+(N3-1)*dim+dim] = hopping_3.transpose().conj() + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 构建有限尺寸的SSH模型哈密顿量 +def hamiltonian_of_finite_size_ssh_model(N, v=0.6, w=1, onsite_1=0, onsite_2=0, period=1): + import numpy as np + hamiltonian = np.zeros((2*N, 2*N)) + for i in range(N): + hamiltonian[i*2+0, i*2+0] = onsite_1 + hamiltonian[i*2+1, i*2+1] = onsite_2 + hamiltonian[i*2+0, i*2+1] = v + hamiltonian[i*2+1, i*2+0] = v + for i in range(N-1): + hamiltonian[i*2+1, (i+1)*2+0] = w + hamiltonian[(i+1)*2+0, i*2+1] = w + if period==1: + hamiltonian[0, 2*N-1] = w + hamiltonian[2*N-1, 0] = w + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 获取Zigzag边的石墨烯条带的元胞间跃迁 +def get_hopping_term_of_graphene_ribbon_along_zigzag_direction(N, eta=0): + import numpy as np + hopping = np.zeros((4*N, 4*N), dtype=complex) + for i0 in range(N): + hopping[4*i0+0, 4*i0+0] = eta + hopping[4*i0+1, 4*i0+1] = eta + hopping[4*i0+2, 4*i0+2] = eta + hopping[4*i0+3, 4*i0+3] = eta + hopping[4*i0+1, 4*i0+0] = 1 + hopping[4*i0+2, 4*i0+3] = 1 + import guan + guan.statistics_of_guan_package() + return hopping + +# 构建有限尺寸的石墨烯哈密顿量（可设置是否为周期边界条件） +def hamiltonian_of_finite_size_system_along_two_directions_for_graphene(N1, N2, period_1=0, period_2=0): + import numpy as np + import guan + on_site = guan.hamiltonian_of_finite_size_system_along_one_direction(4) + hopping_1 = guan.get_hopping_term_of_graphene_ribbon_along_zigzag_direction(1) + hopping_2 = np.zeros((4, 4), dtype=complex) + hopping_2[3, 0] = 1 + hamiltonian = guan.hamiltonian_of_finite_size_system_along_two_directions_for_square_lattice(N1, N2, on_site, hopping_1, hopping_2, period_1, period_2) + guan.statistics_of_guan_package() + return hamiltonian + +# 获取石墨烯有效模型沿着x方向的在位能和跃迁项（其中，动量qy为参数） +def get_onsite_and_hopping_terms_of_2d_effective_graphene_along_one_direction(qy, t=1, staggered_potential=0, eta=0, valley_index=0): + import numpy as np + constant = -np.sqrt(3)/2 + h00 = np.zeros((2, 2), dtype=complex) + h00[0, 0] = staggered_potential + h00[1, 1] = -staggered_potential + h00[0, 1] = -1j*constant*t*np.sin(qy) + h00[1, 0] = 1j*constant*t*np.sin(qy) + h01 = np.zeros((2, 2), dtype=complex) + h01[0, 0] = eta + h01[1, 1] = eta + if valley_index == 0: + h01[0, 1] = constant*t*(-1j/2) + h01[1, 0] = constant*t*(-1j/2) + else: + h01[0, 1] = constant*t*(1j/2) + h01[1, 0] = constant*t*(1j/2) + import guan + guan.statistics_of_guan_package() + return h00, h01 + +# 获取BHZ模型的在位能和跃迁项 +def get_onsite_and_hopping_terms_of_bhz_model(A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01, a=1): + import numpy as np + E_s = C+M-4*(D+B)/(a**2) + E_p = C-M-4*(D-B)/(a**2) + V_ss = (D+B)/(a**2) + V_pp = (D-B)/(a**2) + V_sp = -1j*A/(2*a) + H0 = np.zeros((4, 4), dtype=complex) + H1 = np.zeros((4, 4), dtype=complex) + H2 = np.zeros((4, 4), dtype=complex) + H0[0, 0] = E_s + H0[1, 1] = E_p + H0[2, 2] = E_s + H0[3, 3] = E_p + H1[0, 0] = V_ss + H1[1, 1] = V_pp + H1[2, 2] = V_ss + H1[3, 3] = V_pp + H1[0, 1] = V_sp + H1[1, 0] = -np.conj(V_sp) + H1[2, 3] = np.conj(V_sp) + H1[3, 2] = -V_sp + H2[0, 0] = V_ss + H2[1, 1] = V_pp + H2[2, 2] = V_ss + H2[3, 3] = V_pp + H2[0, 1] = 1j*V_sp + H2[1, 0] = 1j*np.conj(V_sp) + H2[2, 3] = -1j*np.conj(V_sp) + H2[3, 2] = -1j*V_sp + import guan + guan.statistics_of_guan_package() + return H0, H1, H2 + +# 获取半个BHZ模型的在位能和跃迁项（自旋向上） +def get_onsite_and_hopping_terms_of_half_bhz_model_for_spin_up(A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01, a=1): + import numpy as np + E_s = C+M-4*(D+B)/(a**2) + E_p = C-M-4*(D-B)/(a**2) + V_ss = (D+B)/(a**2) + V_pp = (D-B)/(a**2) + V_sp = -1j*A/(2*a) + H0 = np.zeros((2, 2), dtype=complex) + H1 = np.zeros((2, 2), dtype=complex) + H2 = np.zeros((2, 2), dtype=complex) + H0[0, 0] = E_s + H0[1, 1] = E_p + H1[0, 0] = V_ss + H1[1, 1] = V_pp + H1[0, 1] = V_sp + H1[1, 0] = -np.conj(V_sp) + H2[0, 0] = V_ss + H2[1, 1] = V_pp + H2[0, 1] = 1j*V_sp + H2[1, 0] = 1j*np.conj(V_sp) + import guan + guan.statistics_of_guan_package() + return H0, H1, H2 + +# 获取半个BHZ模型的在位能和跃迁项（自旋向下） +def get_onsite_and_hopping_terms_of_half_bhz_model_for_spin_down(A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01, a=1): + import numpy as np + E_s = C+M-4*(D+B)/(a**2) + E_p = C-M-4*(D-B)/(a**2) + V_ss = (D+B)/(a**2) + V_pp = (D-B)/(a**2) + V_sp = -1j*A/(2*a) + H0 = np.zeros((2, 2), dtype=complex) + H1 = np.zeros((2, 2), dtype=complex) + H2 = np.zeros((2, 2), dtype=complex) + H0[0, 0] = E_s + H0[1, 1] = E_p + H1[0, 0] = V_ss + H1[1, 1] = V_pp + H1[0, 1] = np.conj(V_sp) + H1[1, 0] = -V_sp + H2[0, 0] = V_ss + H2[1, 1] = V_pp + H2[0, 1] = -1j*np.conj(V_sp) + H2[1, 0] = -1j*V_sp + import guan + guan.statistics_of_guan_package() + return H0, H1, H2 \ No newline at end of file diff --git a/PyPI/src/guan/Hamiltonian_of_models_in_reciprocal_space.py b/PyPI/src/guan/Hamiltonian_of_models_in_reciprocal_space.py new file mode 100644 index 0000000..ffd5bbd --- /dev/null +++ b/PyPI/src/guan/Hamiltonian_of_models_in_reciprocal_space.py @@ -0,0 +1,315 @@ +# Module: Hamiltonian_of_models_in_reciprocal_space + +# 一维链的哈密顿量 +def hamiltonian_of_simple_chain(k): + import guan + hamiltonian = guan.one_dimensional_fourier_transform(k, unit_cell=0, hopping=1) + guan.statistics_of_guan_package() + return hamiltonian + +# 二维方格子的哈密顿量 +def hamiltonian_of_square_lattice(k1, k2): + import guan + hamiltonian = guan.two_dimensional_fourier_transform_for_square_lattice(k1, k2, unit_cell=0, hopping_1=1, hopping_2=1) + guan.statistics_of_guan_package() + return hamiltonian + +# 准一维方格子条带的哈密顿量 +def hamiltonian_of_square_lattice_in_quasi_one_dimension(k, N=10, period=0): + import numpy as np + import guan + h00 = np.zeros((N, N), dtype=complex) # hopping in a unit cell + h01 = np.zeros((N, N), dtype=complex) # hopping between unit cells + for i in range(N-1): + h00[i, i+1] = 1 + h00[i+1, i] = 1 + if period == 1: + h00[N-1, 0] = 1 + h00[0, N-1] = 1 + for i in range(N): + h01[i, i] = 1 + hamiltonian = guan.one_dimensional_fourier_transform(k, unit_cell=h00, hopping=h01) + guan.statistics_of_guan_package() + return hamiltonian + +# 三维立方格子的哈密顿量 +def hamiltonian_of_cubic_lattice(k1, k2, k3): + import guan + hamiltonian = guan.three_dimensional_fourier_transform_for_cubic_lattice(k1, k2, k3, unit_cell=0, hopping_1=1, hopping_2=1, hopping_3=1) + guan.statistics_of_guan_package() + return hamiltonian + +# SSH模型的哈密顿量 +def hamiltonian_of_ssh_model(k, v=0.6, w=1): + import numpy as np + import cmath + hamiltonian = np.zeros((2, 2), dtype=complex) + hamiltonian[0,1] = v+w*cmath.exp(-1j*k) + hamiltonian[1,0] = v+w*cmath.exp(1j*k) + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 石墨烯的哈密顿量 +def hamiltonian_of_graphene(k1, k2, staggered_potential=0, t=1, a='default'): + import numpy as np + import cmath + import math + if a == 'default': + a = 1/math.sqrt(3) + h0 = np.zeros((2, 2), dtype=complex) # mass term + h1 = np.zeros((2, 2), dtype=complex) # nearest hopping + h0[0, 0] = staggered_potential + h0[1, 1] = -staggered_potential + h1[1, 0] = t*(cmath.exp(1j*k2*a)+cmath.exp(1j*math.sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*math.sqrt(3)/2*k1*a-1j/2*k2*a)) + h1[0, 1] = h1[1, 0].conj() + hamiltonian = h0 + h1 + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 石墨烯有效模型的哈密顿量 +def effective_hamiltonian_of_graphene(qx, qy, t=1, staggered_potential=0, valley_index=0): + import numpy as np + hamiltonian = np.zeros((2, 2), dtype=complex) + hamiltonian[0, 0] = staggered_potential + hamiltonian[1, 1] = -staggered_potential + constant = -np.sqrt(3)/2 + if valley_index == 0: + hamiltonian[0, 1] = constant*t*(qx-1j*qy) + hamiltonian[1, 0] = constant*t*(qx+1j*qy) + else: + hamiltonian[0, 1] = constant*t*(-qx-1j*qy) + hamiltonian[1, 0] = constant*t*(-qx+1j*qy) + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 石墨烯有效模型离散化后的哈密顿量 +def effective_hamiltonian_of_graphene_after_discretization(qx, qy, t=1, staggered_potential=0, valley_index=0): + import numpy as np + hamiltonian = np.zeros((2, 2), dtype=complex) + hamiltonian[0, 0] = staggered_potential + hamiltonian[1, 1] = -staggered_potential + constant = -np.sqrt(3)/2 + if valley_index == 0: + hamiltonian[0, 1] = constant*t*(np.sin(qx)-1j*np.sin(qy)) + hamiltonian[1, 0] = constant*t*(np.sin(qx)+1j*np.sin(qy)) + else: + hamiltonian[0, 1] = constant*t*(-np.sin(qx)-1j*np.sin(qy)) + hamiltonian[1, 0] = constant*t*(-np.sin(qx)+1j*np.sin(qy)) + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 准一维Zigzag边石墨烯条带的哈密顿量 +def hamiltonian_of_graphene_with_zigzag_in_quasi_one_dimension(k, N=10, M=0, t=1, period=0): + import numpy as np + import guan + h00 = np.zeros((4*N, 4*N), dtype=complex) # hopping in a unit cell + h01 = np.zeros((4*N, 4*N), dtype=complex) # hopping between unit cells + for i in range(N): + h00[i*4+0, i*4+0] = M + h00[i*4+1, i*4+1] = -M + h00[i*4+2, i*4+2] = M + h00[i*4+3, i*4+3] = -M + h00[i*4+0, i*4+1] = t + h00[i*4+1, i*4+0] = t + h00[i*4+1, i*4+2] = t + h00[i*4+2, i*4+1] = t + h00[i*4+2, i*4+3] = t + h00[i*4+3, i*4+2] = t + for i in range(N-1): + h00[i*4+3, (i+1)*4+0] = t + h00[(i+1)*4+0, i*4+3] = t + if period == 1: + h00[(N-1)*4+3, 0] = t + h00[0, (N-1)*4+3] = t + for i in range(N): + h01[i*4+1, i*4+0] = t + h01[i*4+2, i*4+3] = t + hamiltonian = guan.one_dimensional_fourier_transform(k, unit_cell=h00, hopping=h01) + guan.statistics_of_guan_package() + return hamiltonian + +# Haldane模型的哈密顿量 +def hamiltonian_of_haldane_model(k1, k2, M=2/3, t1=1, t2=1/3, phi='default', a='default'): + import numpy as np + import cmath + import math + if phi == 'default': + phi=math.pi/4 + if a == 'default': + a=1/math.sqrt(3) + h0 = np.zeros((2, 2), dtype=complex) # mass term + h1 = np.zeros((2, 2), dtype=complex) # nearest hopping + h2 = np.zeros((2, 2), dtype=complex) # next nearest hopping + h0[0, 0] = M + h0[1, 1] = -M + h1[1, 0] = t1*(cmath.exp(1j*k2*a)+cmath.exp(1j*math.sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*math.sqrt(3)/2*k1*a-1j/2*k2*a)) + h1[0, 1] = h1[1, 0].conj() + h2[0, 0] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*math.sqrt(3)*k1*a)+cmath.exp(-1j*math.sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*math.sqrt(3)/2*k1*a-1j*3/2*k2*a)) + h2[1, 1] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*math.sqrt(3)*k1*a)+cmath.exp(-1j*math.sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*math.sqrt(3)/2*k1*a-1j*3/2*k2*a)) + hamiltonian = h0 + h1 + h2 + h2.transpose().conj() + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 准一维Haldane模型条带的哈密顿量 +def hamiltonian_of_haldane_model_in_quasi_one_dimension(k, N=10, M=2/3, t1=1, t2=1/3, phi='default', period=0): + import numpy as np + import cmath + import math + if phi == 'default': + phi=math.pi/4 + h00 = np.zeros((4*N, 4*N), dtype=complex) # hopping in a unit cell + h01 = np.zeros((4*N, 4*N), dtype=complex) # hopping between unit cells + for i in range(N): + h00[i*4+0, i*4+0] = M + h00[i*4+1, i*4+1] = -M + h00[i*4+2, i*4+2] = M + h00[i*4+3, i*4+3] = -M + h00[i*4+0, i*4+1] = t1 + h00[i*4+1, i*4+0] = t1 + h00[i*4+1, i*4+2] = t1 + h00[i*4+2, i*4+1] = t1 + h00[i*4+2, i*4+3] = t1 + h00[i*4+3, i*4+2] = t1 + h00[i*4+0, i*4+2] = t2*cmath.exp(-1j*phi) + h00[i*4+2, i*4+0] = h00[i*4+0, i*4+2].conj() + h00[i*4+1, i*4+3] = t2*cmath.exp(-1j*phi) + h00[i*4+3, i*4+1] = h00[i*4+1, i*4+3].conj() + for i in range(N-1): + h00[i*4+3, (i+1)*4+0] = t1 + h00[(i+1)*4+0, i*4+3] = t1 + h00[i*4+2, (i+1)*4+0] = t2*cmath.exp(1j*phi) + h00[(i+1)*4+0, i*4+2] = h00[i*4+2, (i+1)*4+0].conj() + h00[i*4+3, (i+1)*4+1] = t2*cmath.exp(1j*phi) + h00[(i+1)*4+1, i*4+3] = h00[i*4+3, (i+1)*4+1].conj() + if period == 1: + h00[(N-1)*4+3, 0] = t1 + h00[0, (N-1)*4+3] = t1 + h00[(N-1)*4+2, 0] = t2*cmath.exp(1j*phi) + h00[0, (N-1)*4+2] = h00[(N-1)*4+2, 0].conj() + h00[(N-1)*4+3, 1] = t2*cmath.exp(1j*phi) + h00[1, (N-1)*4+3] = h00[(N-1)*4+3, 1].conj() + for i in range(N): + h01[i*4+1, i*4+0] = t1 + h01[i*4+2, i*4+3] = t1 + h01[i*4+0, i*4+0] = t2*cmath.exp(1j*phi) + h01[i*4+1, i*4+1] = t2*cmath.exp(-1j*phi) + h01[i*4+2, i*4+2] = t2*cmath.exp(1j*phi) + h01[i*4+3, i*4+3] = t2*cmath.exp(-1j*phi) + h01[i*4+1, i*4+3] = t2*cmath.exp(1j*phi) + h01[i*4+2, i*4+0] = t2*cmath.exp(-1j*phi) + if i != 0: + h01[i*4+1, (i-1)*4+3] = t2*cmath.exp(1j*phi) + for i in range(N-1): + h01[i*4+2, (i+1)*4+0] = t2*cmath.exp(-1j*phi) + hamiltonian = h00 + h01*cmath.exp(1j*k) + h01.transpose().conj()*cmath.exp(-1j*k) + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 一个量子反常霍尔效应的哈密顿量 +def hamiltonian_of_one_QAH_model(k1, k2, t1=1, t2=1, t3=0.5, m=-1): + import numpy as np + import math + hamiltonian = np.zeros((2, 2), dtype=complex) + hamiltonian[0, 1] = 2*t1*math.cos(k1)-1j*2*t1*math.cos(k2) + hamiltonian[1, 0] = 2*t1*math.cos(k1)+1j*2*t1*math.cos(k2) + hamiltonian[0, 0] = m+2*t3*math.sin(k1)+2*t3*math.sin(k2)+2*t2*math.cos(k1+k2) + hamiltonian[1, 1] = -(m+2*t3*math.sin(k1)+2*t3*math.sin(k2)+2*t2*math.cos(k1+k2)) + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# BHZ模型的哈密顿量 +def hamiltonian_of_bhz_model(kx, ky, A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01): + import numpy as np + import math + hamiltonian = np.zeros((4, 4), dtype=complex) + varepsilon = C-2*D*(2-math.cos(kx)-math.cos(ky)) + d3 = -2*B*(2-(M/2/B)-math.cos(kx)-math.cos(ky)) + d1_d2 = A*(math.sin(kx)+1j*math.sin(ky)) + hamiltonian[0, 0] = varepsilon+d3 + hamiltonian[1, 1] = varepsilon-d3 + hamiltonian[0, 1] = np.conj(d1_d2) + hamiltonian[1, 0] = d1_d2 + hamiltonian[2, 2] = varepsilon+d3 + hamiltonian[3, 3] = varepsilon-d3 + hamiltonian[2, 3] = -d1_d2 + hamiltonian[3, 2] = -np.conj(d1_d2) + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 半BHZ模型的哈密顿量（自旋向上） +def hamiltonian_of_half_bhz_model_for_spin_up(kx, ky, A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01): + import numpy as np + import math + hamiltonian = np.zeros((2, 2), dtype=complex) + varepsilon = C-2*D*(2-math.cos(kx)-math.cos(ky)) + d3 = -2*B*(2-(M/2/B)-math.cos(kx)-math.cos(ky)) + d1_d2 = A*(math.sin(kx)+1j*math.sin(ky)) + hamiltonian[0, 0] = varepsilon+d3 + hamiltonian[1, 1] = varepsilon-d3 + hamiltonian[0, 1] = np.conj(d1_d2) + hamiltonian[1, 0] = d1_d2 + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# 半BHZ模型的哈密顿量（自旋向下） +def hamiltonian_of_half_bhz_model_for_spin_down(kx, ky, A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01): + import numpy as np + import math + hamiltonian = np.zeros((2, 2), dtype=complex) + varepsilon = C-2*D*(2-math.cos(kx)-math.cos(ky)) + d3 = -2*B*(2-(M/2/B)-math.cos(kx)-math.cos(ky)) + d1_d2 = A*(math.sin(kx)+1j*math.sin(ky)) + hamiltonian[0, 0] = varepsilon+d3 + hamiltonian[1, 1] = varepsilon-d3 + hamiltonian[0, 1] = -d1_d2 + hamiltonian[1, 0] = -np.conj(d1_d2) + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# BBH模型的哈密顿量 +def hamiltonian_of_bbh_model(kx, ky, gamma_x=0.5, gamma_y=0.5, lambda_x=1, lambda_y=1): + import numpy as np + import cmath + # label of atoms in a unit cell + # (2) —— (0) + # | | + # (1) —— (3) + hamiltonian = np.zeros((4, 4), dtype=complex) + hamiltonian[0, 2] = gamma_x+lambda_x*cmath.exp(1j*kx) + hamiltonian[1, 3] = gamma_x+lambda_x*cmath.exp(-1j*kx) + hamiltonian[0, 3] = gamma_y+lambda_y*cmath.exp(1j*ky) + hamiltonian[1, 2] = -gamma_y-lambda_y*cmath.exp(-1j*ky) + hamiltonian[2, 0] = np.conj(hamiltonian[0, 2]) + hamiltonian[3, 1] = np.conj(hamiltonian[1, 3]) + hamiltonian[3, 0] = np.conj(hamiltonian[0, 3]) + hamiltonian[2, 1] = np.conj(hamiltonian[1, 2]) + import guan + guan.statistics_of_guan_package() + return hamiltonian + +# Kagome模型的哈密顿量 +def hamiltonian_of_kagome_lattice(kx, ky, t=1): + import numpy as np + import math + k1_dot_a1 = kx + k2_dot_a2 = kx/2+ky*math.sqrt(3)/2 + k3_dot_a3 = -kx/2+ky*math.sqrt(3)/2 + hamiltonian = np.zeros((3, 3), dtype=complex) + hamiltonian[0, 1] = 2*math.cos(k1_dot_a1) + hamiltonian[0, 2] = 2*math.cos(k2_dot_a2) + hamiltonian[1, 2] = 2*math.cos(k3_dot_a3) + hamiltonian = hamiltonian + hamiltonian.transpose().conj() + hamiltonian = -t*hamiltonian + import guan + guan.statistics_of_guan_package() + return hamiltonian \ No newline at end of file diff --git a/PyPI/src/guan/__init__.py b/PyPI/src/guan/__init__.py index 11e2dee..182b332 100644 --- a/PyPI/src/guan/__init__.py +++ b/PyPI/src/guan/__init__.py @@ -1,4964 +1,20 @@ -# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about (Ji-Huan Guan, 关济寰). The primary location of this package is on website https://py.guanjihuan.com. GitHub link: https://github.com/guanjihuan/py.guanjihuan.com. - -# The current version is guan-0.1.12, updated on November 02, 2023. - -# Installation: pip install --upgrade guan - -# Modules: - -# # Module 1: basic functions -# # Module 2: Fourier transform -# # Module 3: Hamiltonian of finite size systems -# # Module 4: Hamiltonian of models in the reciprocal space -# # Module 5: band structures and wave functions -# # Module 6: Green functions -# # Module 7: density of states -# # Module 8: quantum transport -# # Module 9: topological invariant -# # Module 10: plot figures -# # Module 11: read and write -# # Module 12: file processing -# # Module 13: data processing - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -# Module 1: basic functions - -# 测试 -def test(): - print('\nSuccess in the installation of Guan package!\n') - import guan - guan.statistics_of_guan_package() - -# 泡利矩阵 -def sigma_0(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.eye(2) - -def sigma_x(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.array([[0, 1],[1, 0]]) - -def sigma_y(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.array([[0, -1j],[1j, 0]]) - -def sigma_z(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.array([[1, 0],[0, -1]]) - -# 泡利矩阵的张量积 -def sigma_00(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_0(), guan.sigma_0()) - -def sigma_0x(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_0(), guan.sigma_x()) - -def sigma_0y(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_0(), guan.sigma_y()) - -def sigma_0z(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_0(), guan.sigma_z()) - -def sigma_x0(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_x(), guan.sigma_0()) - -def sigma_xx(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_x(), guan.sigma_x()) - -def sigma_xy(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_x(), guan.sigma_y()) - -def sigma_xz(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_x(), guan.sigma_z()) - -def sigma_y0(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_y(), guan.sigma_0()) - -def sigma_yx(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_y(), guan.sigma_x()) - -def sigma_yy(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_y(), guan.sigma_y()) - -def sigma_yz(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_y(), guan.sigma_z()) - -def sigma_z0(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_z(), guan.sigma_0()) - -def sigma_zx(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_z(), guan.sigma_x()) - -def sigma_zy(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_z(), guan.sigma_y()) - -def sigma_zz(): - import numpy as np - import guan - guan.statistics_of_guan_package() - return np.kron(guan.sigma_z(), guan.sigma_z()) - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -# Module 2: Fourier_transform - -# 通过元胞和跃迁项得到一维的哈密顿量（需要输入k值） -def one_dimensional_fourier_transform(k, unit_cell, hopping): - import numpy as np - import cmath - unit_cell = np.array(unit_cell) - hopping = np.array(hopping) - hamiltonian = unit_cell+hopping*cmath.exp(1j*k)+hopping.transpose().conj()*cmath.exp(-1j*k) - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 通过元胞和跃迁项得到二维方格子的哈密顿量（需要输入k值） -def two_dimensional_fourier_transform_for_square_lattice(k1, k2, unit_cell, hopping_1, hopping_2): - import numpy as np - import cmath - unit_cell = np.array(unit_cell) - hopping_1 = np.array(hopping_1) - hopping_2 = np.array(hopping_2) - hamiltonian = unit_cell+hopping_1*cmath.exp(1j*k1)+hopping_1.transpose().conj()*cmath.exp(-1j*k1)+hopping_2*cmath.exp(1j*k2)+hopping_2.transpose().conj()*cmath.exp(-1j*k2) - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 通过元胞和跃迁项得到三维立方格子的哈密顿量（需要输入k值） -def three_dimensional_fourier_transform_for_cubic_lattice(k1, k2, k3, unit_cell, hopping_1, hopping_2, hopping_3): - import numpy as np - import cmath - unit_cell = np.array(unit_cell) - hopping_1 = np.array(hopping_1) - hopping_2 = np.array(hopping_2) - hopping_3 = np.array(hopping_3) - hamiltonian = unit_cell+hopping_1*cmath.exp(1j*k1)+hopping_1.transpose().conj()*cmath.exp(-1j*k1)+hopping_2*cmath.exp(1j*k2)+hopping_2.transpose().conj()*cmath.exp(-1j*k2)+hopping_3*cmath.exp(1j*k3)+hopping_3.transpose().conj()*cmath.exp(-1j*k3) - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 通过元胞和跃迁项得到一维的哈密顿量（返回的哈密顿量为携带k的函数） -def one_dimensional_fourier_transform_with_k(unit_cell, hopping): - import functools - import guan - hamiltonian_function = functools.partial(guan.one_dimensional_fourier_transform, unit_cell=unit_cell, hopping=hopping) - guan.statistics_of_guan_package() - return hamiltonian_function - -# 通过元胞和跃迁项得到二维方格子的哈密顿量（返回的哈密顿量为携带k的函数） -def two_dimensional_fourier_transform_for_square_lattice_with_k1_k2(unit_cell, hopping_1, hopping_2): - import functools - import guan - hamiltonian_function = functools.partial(guan.two_dimensional_fourier_transform_for_square_lattice, unit_cell=unit_cell, hopping_1=hopping_1, hopping_2=hopping_2) - guan.statistics_of_guan_package() - return hamiltonian_function - -# 通过元胞和跃迁项得到三维立方格子的哈密顿量（返回的哈密顿量为携带k的函数） -def three_dimensional_fourier_transform_for_cubic_lattice_with_k1_k2_k3(unit_cell, hopping_1, hopping_2, hopping_3): - import functools - import guan - hamiltonian_function = functools.partial(guan.three_dimensional_fourier_transform_for_cubic_lattice, unit_cell=unit_cell, hopping_1=hopping_1, hopping_2=hopping_2, hopping_3=hopping_3) - guan.statistics_of_guan_package() - return hamiltonian_function - -# 由实空间格矢得到倒空间格矢（一维） -def calculate_one_dimensional_reciprocal_lattice_vector(a1): - import numpy as np - b1 = 2*np.pi/a1 - import guan - guan.statistics_of_guan_package() - return b1 - -# 由实空间格矢得到倒空间格矢（二维） -def calculate_two_dimensional_reciprocal_lattice_vectors(a1, a2): - import numpy as np - a1 = np.array(a1) - a2 = np.array(a2) - a1 = np.append(a1, 0) - a2 = np.append(a2, 0) - a3 = np.array([0, 0, 1]) - b1 = 2*np.pi*np.cross(a2, a3)/np.dot(a1, np.cross(a2, a3)) - b2 = 2*np.pi*np.cross(a3, a1)/np.dot(a1, np.cross(a2, a3)) - b1 = np.delete(b1, 2) - b2 = np.delete(b2, 2) - import guan - guan.statistics_of_guan_package() - return b1, b2 - -# 由实空间格矢得到倒空间格矢（三维） -def calculate_three_dimensional_reciprocal_lattice_vectors(a1, a2, a3): - import numpy as np - a1 = np.array(a1) - a2 = np.array(a2) - a3 = np.array(a3) - b1 = 2*np.pi*np.cross(a2, a3)/np.dot(a1, np.cross(a2, a3)) - b2 = 2*np.pi*np.cross(a3, a1)/np.dot(a1, np.cross(a2, a3)) - b3 = 2*np.pi*np.cross(a1, a2)/np.dot(a1, np.cross(a2, a3)) - import guan - guan.statistics_of_guan_package() - return b1, b2, b3 - -# 由实空间格矢得到倒空间格矢（一维），这里为符号运算 -def calculate_one_dimensional_reciprocal_lattice_vector_with_sympy(a1): - import sympy - b1 = 2*sympy.pi/a1 - import guan - guan.statistics_of_guan_package() - return b1 - -# 由实空间格矢得到倒空间格矢（二维），这里为符号运算 -def calculate_two_dimensional_reciprocal_lattice_vectors_with_sympy(a1, a2): - import sympy - a1 = sympy.Matrix(1, 3, [a1[0], a1[1], 0]) - a2 = sympy.Matrix(1, 3, [a2[0], a2[1], 0]) - a3 = sympy.Matrix(1, 3, [0, 0, 1]) - cross_a2_a3 = a2.cross(a3) - cross_a3_a1 = a3.cross(a1) - b1 = 2*sympy.pi*cross_a2_a3/a1.dot(cross_a2_a3) - b2 = 2*sympy.pi*cross_a3_a1/a1.dot(cross_a2_a3) - b1 = sympy.Matrix(1, 2, [b1[0], b1[1]]) - b2 = sympy.Matrix(1, 2, [b2[0], b2[1]]) - import guan - guan.statistics_of_guan_package() - return b1, b2 - -# 由实空间格矢得到倒空间格矢（三维），这里为符号运算 -def calculate_three_dimensional_reciprocal_lattice_vectors_with_sympy(a1, a2, a3): - import sympy - cross_a2_a3 = a2.cross(a3) - cross_a3_a1 = a3.cross(a1) - cross_a1_a2 = a1.cross(a2) - b1 = 2*sympy.pi*cross_a2_a3/a1.dot(cross_a2_a3) - b2 = 2*sympy.pi*cross_a3_a1/a1.dot(cross_a2_a3) - b3 = 2*sympy.pi*cross_a1_a2/a1.dot(cross_a2_a3) - import guan - guan.statistics_of_guan_package() - return b1, b2, b3 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -# Module 3: Hamiltonian of finite size systems - -# 构建一维的有限尺寸体系哈密顿量（可设置是否为周期边界条件） -def hamiltonian_of_finite_size_system_along_one_direction(N, on_site=0, hopping=1, period=0): - import numpy as np - on_site = np.array(on_site) - hopping = np.array(hopping) - if on_site.shape==(): - dim = 1 - else: - dim = on_site.shape[0] - hamiltonian = np.zeros((N*dim, N*dim), dtype=complex) - for i0 in range(N): - hamiltonian[i0*dim+0:i0*dim+dim, i0*dim+0:i0*dim+dim] = on_site - for i0 in range(N-1): - hamiltonian[i0*dim+0:i0*dim+dim, (i0+1)*dim+0:(i0+1)*dim+dim] = hopping - hamiltonian[(i0+1)*dim+0:(i0+1)*dim+dim, i0*dim+0:i0*dim+dim] = hopping.transpose().conj() - if period == 1: - hamiltonian[(N-1)*dim+0:(N-1)*dim+dim, 0:dim] = hopping - hamiltonian[0:dim, (N-1)*dim+0:(N-1)*dim+dim] = hopping.transpose().conj() - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 构建二维的方格子有限尺寸体系哈密顿量（可设置是否为周期边界条件） -def hamiltonian_of_finite_size_system_along_two_directions_for_square_lattice(N1, N2, on_site=0, hopping_1=1, hopping_2=1, period_1=0, period_2=0): - import numpy as np - on_site = np.array(on_site) - hopping_1 = np.array(hopping_1) - hopping_2 = np.array(hopping_2) - if on_site.shape==(): - dim = 1 - else: - dim = on_site.shape[0] - hamiltonian = np.zeros((N1*N2*dim, N1*N2*dim), dtype=complex) - for i1 in range(N1): - for i2 in range(N2): - hamiltonian[i1*N2*dim+i2*dim+0:i1*N2*dim+i2*dim+dim, i1*N2*dim+i2*dim+0:i1*N2*dim+i2*dim+dim] = on_site - for i1 in range(N1-1): - for i2 in range(N2): - hamiltonian[i1*N2*dim+i2*dim+0:i1*N2*dim+i2*dim+dim, (i1+1)*N2*dim+i2*dim+0:(i1+1)*N2*dim+i2*dim+dim] = hopping_1 - hamiltonian[(i1+1)*N2*dim+i2*dim+0:(i1+1)*N2*dim+i2*dim+dim, i1*N2*dim+i2*dim+0:i1*N2*dim+i2*dim+dim] = hopping_1.transpose().conj() - for i1 in range(N1): - for i2 in range(N2-1): - hamiltonian[i1*N2*dim+i2*dim+0:i1*N2*dim+i2*dim+dim, i1*N2*dim+(i2+1)*dim+0:i1*N2*dim+(i2+1)*dim+dim] = hopping_2 - hamiltonian[i1*N2*dim+(i2+1)*dim+0:i1*N2*dim+(i2+1)*dim+dim, i1*N2*dim+i2*dim+0:i1*N2*dim+i2*dim+dim] = hopping_2.transpose().conj() - if period_1 == 1: - for i2 in range(N2): - hamiltonian[(N1-1)*N2*dim+i2*dim+0:(N1-1)*N2*dim+i2*dim+dim, i2*dim+0:i2*dim+dim] = hopping_1 - hamiltonian[i2*dim+0:i2*dim+dim, (N1-1)*N2*dim+i2*dim+0:(N1-1)*N2*dim+i2*dim+dim] = hopping_1.transpose().conj() - if period_2 == 1: - for i1 in range(N1): - hamiltonian[i1*N2*dim+(N2-1)*dim+0:i1*N2*dim+(N2-1)*dim+dim, i1*N2*dim+0:i1*N2*dim+dim] = hopping_2 - hamiltonian[i1*N2*dim+0:i1*N2*dim+dim, i1*N2*dim+(N2-1)*dim+0:i1*N2*dim+(N2-1)*dim+dim] = hopping_2.transpose().conj() - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 构建三维的立方格子有限尺寸体系哈密顿量（可设置是否为周期边界条件） -def hamiltonian_of_finite_size_system_along_three_directions_for_cubic_lattice(N1, N2, N3, on_site=0, hopping_1=1, hopping_2=1, hopping_3=1, period_1=0, period_2=0, period_3=0): - import numpy as np - on_site = np.array(on_site) - hopping_1 = np.array(hopping_1) - hopping_2 = np.array(hopping_2) - hopping_3 = np.array(hopping_3) - if on_site.shape==(): - dim = 1 - else: - dim = on_site.shape[0] - hamiltonian = np.zeros((N1*N2*N3*dim, N1*N2*N3*dim), dtype=complex) - for i1 in range(N1): - for i2 in range(N2): - for i3 in range(N3): - hamiltonian[i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim, i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim] = on_site - for i1 in range(N1-1): - for i2 in range(N2): - for i3 in range(N3): - hamiltonian[i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim, (i1+1)*N2*N3*dim+i2*N3*dim+i3*dim+0:(i1+1)*N2*N3*dim+i2*N3*dim+i3*dim+dim] = hopping_1 - hamiltonian[(i1+1)*N2*N3*dim+i2*N3*dim+i3*dim+0:(i1+1)*N2*N3*dim+i2*N3*dim+i3*dim+dim, i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim] = hopping_1.transpose().conj() - for i1 in range(N1): - for i2 in range(N2-1): - for i3 in range(N3): - hamiltonian[i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim, i1*N2*N3*dim+(i2+1)*N3*dim+i3*dim+0:i1*N2*N3*dim+(i2+1)*N3*dim+i3*dim+dim] = hopping_2 - hamiltonian[i1*N2*N3*dim+(i2+1)*N3*dim+i3*dim+0:i1*N2*N3*dim+(i2+1)*N3*dim+i3*dim+dim, i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim] = hopping_2.transpose().conj() - for i1 in range(N1): - for i2 in range(N2): - for i3 in range(N3-1): - hamiltonian[i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim, i1*N2*N3*dim+i2*N3*dim+(i3+1)*dim+0:i1*N2*N3*dim+i2*N3*dim+(i3+1)*dim+dim] = hopping_3 - hamiltonian[i1*N2*N3*dim+i2*N3*dim+(i3+1)*dim+0:i1*N2*N3*dim+i2*N3*dim+(i3+1)*dim+dim, i1*N2*N3*dim+i2*N3*dim+i3*dim+0:i1*N2*N3*dim+i2*N3*dim+i3*dim+dim] = hopping_3.transpose().conj() - if period_1 == 1: - for i2 in range(N2): - for i3 in range(N3): - hamiltonian[(N1-1)*N2*N3*dim+i2*N3*dim+i3*dim+0:(N1-1)*N2*N3*dim+i2*N3*dim+i3*dim+dim, i2*N3*dim+i3*dim+0:i2*N3*dim+i3*dim+dim] = hopping_1 - hamiltonian[i2*N3*dim+i3*dim+0:i2*N3*dim+i3*dim+dim, (N1-1)*N2*N3*dim+i2*N3*dim+i3*dim+0:(N1-1)*N2*N3*dim+i2*N3*dim+i3*dim+dim] = hopping_1.transpose().conj() - if period_2 == 1: - for i1 in range(N1): - for i3 in range(N3): - hamiltonian[i1*N2*N3*dim+(N2-1)*N3*dim+i3*dim+0:i1*N2*N3*dim+(N2-1)*N3*dim+i3*dim+dim, i1*N2*N3*dim+i3*dim+0:i1*N2*N3*dim+i3*dim+dim] = hopping_2 - hamiltonian[i1*N2*N3*dim+i3*dim+0:i1*N2*N3*dim+i3*dim+dim, i1*N2*N3*dim+(N2-1)*N3*dim+i3*dim+0:i1*N2*N3*dim+(N2-1)*N3*dim+i3*dim+dim] = hopping_2.transpose().conj() - if period_3 == 1: - for i1 in range(N1): - for i2 in range(N2): - hamiltonian[i1*N2*N3*dim+i2*N3*dim+(N3-1)*dim+0:i1*N2*N3*dim+i2*N3*dim+(N3-1)*dim+dim, i1*N2*N3*dim+i2*N3*dim+0:i1*N2*N3*dim+i2*N3*dim+dim] = hopping_3 - hamiltonian[i1*N2*N3*dim+i2*N3*dim+0:i1*N2*N3*dim+i2*N3*dim+dim, i1*N2*N3*dim+i2*N3*dim+(N3-1)*dim+0:i1*N2*N3*dim+i2*N3*dim+(N3-1)*dim+dim] = hopping_3.transpose().conj() - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 构建有限尺寸的SSH模型哈密顿量 -def hamiltonian_of_finite_size_ssh_model(N, v=0.6, w=1, onsite_1=0, onsite_2=0, period=1): - import numpy as np - hamiltonian = np.zeros((2*N, 2*N)) - for i in range(N): - hamiltonian[i*2+0, i*2+0] = onsite_1 - hamiltonian[i*2+1, i*2+1] = onsite_2 - hamiltonian[i*2+0, i*2+1] = v - hamiltonian[i*2+1, i*2+0] = v - for i in range(N-1): - hamiltonian[i*2+1, (i+1)*2+0] = w - hamiltonian[(i+1)*2+0, i*2+1] = w - if period==1: - hamiltonian[0, 2*N-1] = w - hamiltonian[2*N-1, 0] = w - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 获取Zigzag边的石墨烯条带的元胞间跃迁 -def get_hopping_term_of_graphene_ribbon_along_zigzag_direction(N, eta=0): - import numpy as np - hopping = np.zeros((4*N, 4*N), dtype=complex) - for i0 in range(N): - hopping[4*i0+0, 4*i0+0] = eta - hopping[4*i0+1, 4*i0+1] = eta - hopping[4*i0+2, 4*i0+2] = eta - hopping[4*i0+3, 4*i0+3] = eta - hopping[4*i0+1, 4*i0+0] = 1 - hopping[4*i0+2, 4*i0+3] = 1 - import guan - guan.statistics_of_guan_package() - return hopping - -# 构建有限尺寸的石墨烯哈密顿量（可设置是否为周期边界条件） -def hamiltonian_of_finite_size_system_along_two_directions_for_graphene(N1, N2, period_1=0, period_2=0): - import numpy as np - import guan - on_site = guan.hamiltonian_of_finite_size_system_along_one_direction(4) - hopping_1 = guan.get_hopping_term_of_graphene_ribbon_along_zigzag_direction(1) - hopping_2 = np.zeros((4, 4), dtype=complex) - hopping_2[3, 0] = 1 - hamiltonian = guan.hamiltonian_of_finite_size_system_along_two_directions_for_square_lattice(N1, N2, on_site, hopping_1, hopping_2, period_1, period_2) - guan.statistics_of_guan_package() - return hamiltonian - -# 获取石墨烯有效模型沿着x方向的在位能和跃迁项（其中，动量qy为参数） -def get_onsite_and_hopping_terms_of_2d_effective_graphene_along_one_direction(qy, t=1, staggered_potential=0, eta=0, valley_index=0): - import numpy as np - constant = -np.sqrt(3)/2 - h00 = np.zeros((2, 2), dtype=complex) - h00[0, 0] = staggered_potential - h00[1, 1] = -staggered_potential - h00[0, 1] = -1j*constant*t*np.sin(qy) - h00[1, 0] = 1j*constant*t*np.sin(qy) - h01 = np.zeros((2, 2), dtype=complex) - h01[0, 0] = eta - h01[1, 1] = eta - if valley_index == 0: - h01[0, 1] = constant*t*(-1j/2) - h01[1, 0] = constant*t*(-1j/2) - else: - h01[0, 1] = constant*t*(1j/2) - h01[1, 0] = constant*t*(1j/2) - import guan - guan.statistics_of_guan_package() - return h00, h01 - -# 获取BHZ模型的在位能和跃迁项 -def get_onsite_and_hopping_terms_of_bhz_model(A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01, a=1): - import numpy as np - E_s = C+M-4*(D+B)/(a**2) - E_p = C-M-4*(D-B)/(a**2) - V_ss = (D+B)/(a**2) - V_pp = (D-B)/(a**2) - V_sp = -1j*A/(2*a) - H0 = np.zeros((4, 4), dtype=complex) - H1 = np.zeros((4, 4), dtype=complex) - H2 = np.zeros((4, 4), dtype=complex) - H0[0, 0] = E_s - H0[1, 1] = E_p - H0[2, 2] = E_s - H0[3, 3] = E_p - H1[0, 0] = V_ss - H1[1, 1] = V_pp - H1[2, 2] = V_ss - H1[3, 3] = V_pp - H1[0, 1] = V_sp - H1[1, 0] = -np.conj(V_sp) - H1[2, 3] = np.conj(V_sp) - H1[3, 2] = -V_sp - H2[0, 0] = V_ss - H2[1, 1] = V_pp - H2[2, 2] = V_ss - H2[3, 3] = V_pp - H2[0, 1] = 1j*V_sp - H2[1, 0] = 1j*np.conj(V_sp) - H2[2, 3] = -1j*np.conj(V_sp) - H2[3, 2] = -1j*V_sp - import guan - guan.statistics_of_guan_package() - return H0, H1, H2 - -# 获取半个BHZ模型的在位能和跃迁项（自旋向上） -def get_onsite_and_hopping_terms_of_half_bhz_model_for_spin_up(A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01, a=1): - import numpy as np - E_s = C+M-4*(D+B)/(a**2) - E_p = C-M-4*(D-B)/(a**2) - V_ss = (D+B)/(a**2) - V_pp = (D-B)/(a**2) - V_sp = -1j*A/(2*a) - H0 = np.zeros((2, 2), dtype=complex) - H1 = np.zeros((2, 2), dtype=complex) - H2 = np.zeros((2, 2), dtype=complex) - H0[0, 0] = E_s - H0[1, 1] = E_p - H1[0, 0] = V_ss - H1[1, 1] = V_pp - H1[0, 1] = V_sp - H1[1, 0] = -np.conj(V_sp) - H2[0, 0] = V_ss - H2[1, 1] = V_pp - H2[0, 1] = 1j*V_sp - H2[1, 0] = 1j*np.conj(V_sp) - import guan - guan.statistics_of_guan_package() - return H0, H1, H2 - -# 获取半个BHZ模型的在位能和跃迁项（自旋向下） -def get_onsite_and_hopping_terms_of_half_bhz_model_for_spin_down(A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01, a=1): - import numpy as np - E_s = C+M-4*(D+B)/(a**2) - E_p = C-M-4*(D-B)/(a**2) - V_ss = (D+B)/(a**2) - V_pp = (D-B)/(a**2) - V_sp = -1j*A/(2*a) - H0 = np.zeros((2, 2), dtype=complex) - H1 = np.zeros((2, 2), dtype=complex) - H2 = np.zeros((2, 2), dtype=complex) - H0[0, 0] = E_s - H0[1, 1] = E_p - H1[0, 0] = V_ss - H1[1, 1] = V_pp - H1[0, 1] = np.conj(V_sp) - H1[1, 0] = -V_sp - H2[0, 0] = V_ss - H2[1, 1] = V_pp - H2[0, 1] = -1j*np.conj(V_sp) - H2[1, 0] = -1j*V_sp - import guan - guan.statistics_of_guan_package() - return H0, H1, H2 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -# Module 4: Hamiltonian of models in the reciprocal space - -# 一维链的哈密顿量 -def hamiltonian_of_simple_chain(k): - import guan - hamiltonian = guan.one_dimensional_fourier_transform(k, unit_cell=0, hopping=1) - guan.statistics_of_guan_package() - return hamiltonian - -# 二维方格子的哈密顿量 -def hamiltonian_of_square_lattice(k1, k2): - import guan - hamiltonian = guan.two_dimensional_fourier_transform_for_square_lattice(k1, k2, unit_cell=0, hopping_1=1, hopping_2=1) - guan.statistics_of_guan_package() - return hamiltonian - -# 准一维方格子条带的哈密顿量 -def hamiltonian_of_square_lattice_in_quasi_one_dimension(k, N=10, period=0): - import numpy as np - import guan - h00 = np.zeros((N, N), dtype=complex) # hopping in a unit cell - h01 = np.zeros((N, N), dtype=complex) # hopping between unit cells - for i in range(N-1): - h00[i, i+1] = 1 - h00[i+1, i] = 1 - if period == 1: - h00[N-1, 0] = 1 - h00[0, N-1] = 1 - for i in range(N): - h01[i, i] = 1 - hamiltonian = guan.one_dimensional_fourier_transform(k, unit_cell=h00, hopping=h01) - guan.statistics_of_guan_package() - return hamiltonian - -# 三维立方格子的哈密顿量 -def hamiltonian_of_cubic_lattice(k1, k2, k3): - import guan - hamiltonian = guan.three_dimensional_fourier_transform_for_cubic_lattice(k1, k2, k3, unit_cell=0, hopping_1=1, hopping_2=1, hopping_3=1) - guan.statistics_of_guan_package() - return hamiltonian - -# SSH模型的哈密顿量 -def hamiltonian_of_ssh_model(k, v=0.6, w=1): - import numpy as np - import cmath - hamiltonian = np.zeros((2, 2), dtype=complex) - hamiltonian[0,1] = v+w*cmath.exp(-1j*k) - hamiltonian[1,0] = v+w*cmath.exp(1j*k) - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 石墨烯的哈密顿量 -def hamiltonian_of_graphene(k1, k2, staggered_potential=0, t=1, a='default'): - import numpy as np - import cmath - import math - if a == 'default': - a = 1/math.sqrt(3) - h0 = np.zeros((2, 2), dtype=complex) # mass term - h1 = np.zeros((2, 2), dtype=complex) # nearest hopping - h0[0, 0] = staggered_potential - h0[1, 1] = -staggered_potential - h1[1, 0] = t*(cmath.exp(1j*k2*a)+cmath.exp(1j*math.sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*math.sqrt(3)/2*k1*a-1j/2*k2*a)) - h1[0, 1] = h1[1, 0].conj() - hamiltonian = h0 + h1 - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 石墨烯有效模型的哈密顿量 -def effective_hamiltonian_of_graphene(qx, qy, t=1, staggered_potential=0, valley_index=0): - import numpy as np - hamiltonian = np.zeros((2, 2), dtype=complex) - hamiltonian[0, 0] = staggered_potential - hamiltonian[1, 1] = -staggered_potential - constant = -np.sqrt(3)/2 - if valley_index == 0: - hamiltonian[0, 1] = constant*t*(qx-1j*qy) - hamiltonian[1, 0] = constant*t*(qx+1j*qy) - else: - hamiltonian[0, 1] = constant*t*(-qx-1j*qy) - hamiltonian[1, 0] = constant*t*(-qx+1j*qy) - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 石墨烯有效模型离散化后的哈密顿量 -def effective_hamiltonian_of_graphene_after_discretization(qx, qy, t=1, staggered_potential=0, valley_index=0): - import numpy as np - hamiltonian = np.zeros((2, 2), dtype=complex) - hamiltonian[0, 0] = staggered_potential - hamiltonian[1, 1] = -staggered_potential - constant = -np.sqrt(3)/2 - if valley_index == 0: - hamiltonian[0, 1] = constant*t*(np.sin(qx)-1j*np.sin(qy)) - hamiltonian[1, 0] = constant*t*(np.sin(qx)+1j*np.sin(qy)) - else: - hamiltonian[0, 1] = constant*t*(-np.sin(qx)-1j*np.sin(qy)) - hamiltonian[1, 0] = constant*t*(-np.sin(qx)+1j*np.sin(qy)) - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 准一维Zigzag边石墨烯条带的哈密顿量 -def hamiltonian_of_graphene_with_zigzag_in_quasi_one_dimension(k, N=10, M=0, t=1, period=0): - import numpy as np - import guan - h00 = np.zeros((4*N, 4*N), dtype=complex) # hopping in a unit cell - h01 = np.zeros((4*N, 4*N), dtype=complex) # hopping between unit cells - for i in range(N): - h00[i*4+0, i*4+0] = M - h00[i*4+1, i*4+1] = -M - h00[i*4+2, i*4+2] = M - h00[i*4+3, i*4+3] = -M - h00[i*4+0, i*4+1] = t - h00[i*4+1, i*4+0] = t - h00[i*4+1, i*4+2] = t - h00[i*4+2, i*4+1] = t - h00[i*4+2, i*4+3] = t - h00[i*4+3, i*4+2] = t - for i in range(N-1): - h00[i*4+3, (i+1)*4+0] = t - h00[(i+1)*4+0, i*4+3] = t - if period == 1: - h00[(N-1)*4+3, 0] = t - h00[0, (N-1)*4+3] = t - for i in range(N): - h01[i*4+1, i*4+0] = t - h01[i*4+2, i*4+3] = t - hamiltonian = guan.one_dimensional_fourier_transform(k, unit_cell=h00, hopping=h01) - guan.statistics_of_guan_package() - return hamiltonian - -# Haldane模型的哈密顿量 -def hamiltonian_of_haldane_model(k1, k2, M=2/3, t1=1, t2=1/3, phi='default', a='default'): - import numpy as np - import cmath - import math - if phi == 'default': - phi=math.pi/4 - if a == 'default': - a=1/math.sqrt(3) - h0 = np.zeros((2, 2), dtype=complex) # mass term - h1 = np.zeros((2, 2), dtype=complex) # nearest hopping - h2 = np.zeros((2, 2), dtype=complex) # next nearest hopping - h0[0, 0] = M - h0[1, 1] = -M - h1[1, 0] = t1*(cmath.exp(1j*k2*a)+cmath.exp(1j*math.sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*math.sqrt(3)/2*k1*a-1j/2*k2*a)) - h1[0, 1] = h1[1, 0].conj() - h2[0, 0] = t2*cmath.exp(-1j*phi)*(cmath.exp(1j*math.sqrt(3)*k1*a)+cmath.exp(-1j*math.sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*math.sqrt(3)/2*k1*a-1j*3/2*k2*a)) - h2[1, 1] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*math.sqrt(3)*k1*a)+cmath.exp(-1j*math.sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*math.sqrt(3)/2*k1*a-1j*3/2*k2*a)) - hamiltonian = h0 + h1 + h2 + h2.transpose().conj() - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 准一维Haldane模型条带的哈密顿量 -def hamiltonian_of_haldane_model_in_quasi_one_dimension(k, N=10, M=2/3, t1=1, t2=1/3, phi='default', period=0): - import numpy as np - import cmath - import math - if phi == 'default': - phi=math.pi/4 - h00 = np.zeros((4*N, 4*N), dtype=complex) # hopping in a unit cell - h01 = np.zeros((4*N, 4*N), dtype=complex) # hopping between unit cells - for i in range(N): - h00[i*4+0, i*4+0] = M - h00[i*4+1, i*4+1] = -M - h00[i*4+2, i*4+2] = M - h00[i*4+3, i*4+3] = -M - h00[i*4+0, i*4+1] = t1 - h00[i*4+1, i*4+0] = t1 - h00[i*4+1, i*4+2] = t1 - h00[i*4+2, i*4+1] = t1 - h00[i*4+2, i*4+3] = t1 - h00[i*4+3, i*4+2] = t1 - h00[i*4+0, i*4+2] = t2*cmath.exp(-1j*phi) - h00[i*4+2, i*4+0] = h00[i*4+0, i*4+2].conj() - h00[i*4+1, i*4+3] = t2*cmath.exp(-1j*phi) - h00[i*4+3, i*4+1] = h00[i*4+1, i*4+3].conj() - for i in range(N-1): - h00[i*4+3, (i+1)*4+0] = t1 - h00[(i+1)*4+0, i*4+3] = t1 - h00[i*4+2, (i+1)*4+0] = t2*cmath.exp(1j*phi) - h00[(i+1)*4+0, i*4+2] = h00[i*4+2, (i+1)*4+0].conj() - h00[i*4+3, (i+1)*4+1] = t2*cmath.exp(1j*phi) - h00[(i+1)*4+1, i*4+3] = h00[i*4+3, (i+1)*4+1].conj() - if period == 1: - h00[(N-1)*4+3, 0] = t1 - h00[0, (N-1)*4+3] = t1 - h00[(N-1)*4+2, 0] = t2*cmath.exp(1j*phi) - h00[0, (N-1)*4+2] = h00[(N-1)*4+2, 0].conj() - h00[(N-1)*4+3, 1] = t2*cmath.exp(1j*phi) - h00[1, (N-1)*4+3] = h00[(N-1)*4+3, 1].conj() - for i in range(N): - h01[i*4+1, i*4+0] = t1 - h01[i*4+2, i*4+3] = t1 - h01[i*4+0, i*4+0] = t2*cmath.exp(1j*phi) - h01[i*4+1, i*4+1] = t2*cmath.exp(-1j*phi) - h01[i*4+2, i*4+2] = t2*cmath.exp(1j*phi) - h01[i*4+3, i*4+3] = t2*cmath.exp(-1j*phi) - h01[i*4+1, i*4+3] = t2*cmath.exp(1j*phi) - h01[i*4+2, i*4+0] = t2*cmath.exp(-1j*phi) - if i != 0: - h01[i*4+1, (i-1)*4+3] = t2*cmath.exp(1j*phi) - for i in range(N-1): - h01[i*4+2, (i+1)*4+0] = t2*cmath.exp(-1j*phi) - hamiltonian = h00 + h01*cmath.exp(1j*k) + h01.transpose().conj()*cmath.exp(-1j*k) - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 一个量子反常霍尔效应的哈密顿量 -def hamiltonian_of_one_QAH_model(k1, k2, t1=1, t2=1, t3=0.5, m=-1): - import numpy as np - import math - hamiltonian = np.zeros((2, 2), dtype=complex) - hamiltonian[0, 1] = 2*t1*math.cos(k1)-1j*2*t1*math.cos(k2) - hamiltonian[1, 0] = 2*t1*math.cos(k1)+1j*2*t1*math.cos(k2) - hamiltonian[0, 0] = m+2*t3*math.sin(k1)+2*t3*math.sin(k2)+2*t2*math.cos(k1+k2) - hamiltonian[1, 1] = -(m+2*t3*math.sin(k1)+2*t3*math.sin(k2)+2*t2*math.cos(k1+k2)) - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# BHZ模型的哈密顿量 -def hamiltonian_of_bhz_model(kx, ky, A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01): - import numpy as np - import math - hamiltonian = np.zeros((4, 4), dtype=complex) - varepsilon = C-2*D*(2-math.cos(kx)-math.cos(ky)) - d3 = -2*B*(2-(M/2/B)-math.cos(kx)-math.cos(ky)) - d1_d2 = A*(math.sin(kx)+1j*math.sin(ky)) - hamiltonian[0, 0] = varepsilon+d3 - hamiltonian[1, 1] = varepsilon-d3 - hamiltonian[0, 1] = np.conj(d1_d2) - hamiltonian[1, 0] = d1_d2 - hamiltonian[2, 2] = varepsilon+d3 - hamiltonian[3, 3] = varepsilon-d3 - hamiltonian[2, 3] = -d1_d2 - hamiltonian[3, 2] = -np.conj(d1_d2) - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 半BHZ模型的哈密顿量（自旋向上） -def hamiltonian_of_half_bhz_model_for_spin_up(kx, ky, A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01): - import numpy as np - import math - hamiltonian = np.zeros((2, 2), dtype=complex) - varepsilon = C-2*D*(2-math.cos(kx)-math.cos(ky)) - d3 = -2*B*(2-(M/2/B)-math.cos(kx)-math.cos(ky)) - d1_d2 = A*(math.sin(kx)+1j*math.sin(ky)) - hamiltonian[0, 0] = varepsilon+d3 - hamiltonian[1, 1] = varepsilon-d3 - hamiltonian[0, 1] = np.conj(d1_d2) - hamiltonian[1, 0] = d1_d2 - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# 半BHZ模型的哈密顿量（自旋向下） -def hamiltonian_of_half_bhz_model_for_spin_down(kx, ky, A=0.3645/5, B=-0.686/25, C=0, D=-0.512/25, M=-0.01): - import numpy as np - import math - hamiltonian = np.zeros((2, 2), dtype=complex) - varepsilon = C-2*D*(2-math.cos(kx)-math.cos(ky)) - d3 = -2*B*(2-(M/2/B)-math.cos(kx)-math.cos(ky)) - d1_d2 = A*(math.sin(kx)+1j*math.sin(ky)) - hamiltonian[0, 0] = varepsilon+d3 - hamiltonian[1, 1] = varepsilon-d3 - hamiltonian[0, 1] = -d1_d2 - hamiltonian[1, 0] = -np.conj(d1_d2) - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# BBH模型的哈密顿量 -def hamiltonian_of_bbh_model(kx, ky, gamma_x=0.5, gamma_y=0.5, lambda_x=1, lambda_y=1): - import numpy as np - import cmath - # label of atoms in a unit cell - # (2) —— (0) - # | | - # (1) —— (3) - hamiltonian = np.zeros((4, 4), dtype=complex) - hamiltonian[0, 2] = gamma_x+lambda_x*cmath.exp(1j*kx) - hamiltonian[1, 3] = gamma_x+lambda_x*cmath.exp(-1j*kx) - hamiltonian[0, 3] = gamma_y+lambda_y*cmath.exp(1j*ky) - hamiltonian[1, 2] = -gamma_y-lambda_y*cmath.exp(-1j*ky) - hamiltonian[2, 0] = np.conj(hamiltonian[0, 2]) - hamiltonian[3, 1] = np.conj(hamiltonian[1, 3]) - hamiltonian[3, 0] = np.conj(hamiltonian[0, 3]) - hamiltonian[2, 1] = np.conj(hamiltonian[1, 2]) - import guan - guan.statistics_of_guan_package() - return hamiltonian - -# Kagome模型的哈密顿量 -def hamiltonian_of_kagome_lattice(kx, ky, t=1): - import numpy as np - import math - k1_dot_a1 = kx - k2_dot_a2 = kx/2+ky*math.sqrt(3)/2 - k3_dot_a3 = -kx/2+ky*math.sqrt(3)/2 - hamiltonian = np.zeros((3, 3), dtype=complex) - hamiltonian[0, 1] = 2*math.cos(k1_dot_a1) - hamiltonian[0, 2] = 2*math.cos(k2_dot_a2) - hamiltonian[1, 2] = 2*math.cos(k3_dot_a3) - hamiltonian = hamiltonian + hamiltonian.transpose().conj() - hamiltonian = -t*hamiltonian - import guan - guan.statistics_of_guan_package() - return hamiltonian - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -# Module 5: band structures and wave functions - -# 计算哈密顿量的本征值 -def calculate_eigenvalue(hamiltonian): - import numpy as np - if np.array(hamiltonian).shape==(): - eigenvalue = np.real(hamiltonian) - else: - eigenvalue, eigenvector = np.linalg.eigh(hamiltonian) - import guan - guan.statistics_of_guan_package() - return eigenvalue - -# 输入哈密顿量函数（带一组参数），计算一组参数下的本征值，返回本征值向量组 -def calculate_eigenvalue_with_one_parameter(x_array, hamiltonian_function, print_show=0): - import numpy as np - dim_x = np.array(x_array).shape[0] - i0 = 0 - if np.array(hamiltonian_function(0)).shape==(): - eigenvalue_array = np.zeros((dim_x, 1)) - for x0 in x_array: - hamiltonian = hamiltonian_function(x0) - eigenvalue_array[i0, 0] = np.real(hamiltonian) - i0 += 1 - else: - dim = np.array(hamiltonian_function(0)).shape[0] - eigenvalue_array = np.zeros((dim_x, dim)) - for x0 in x_array: - if print_show==1: - print(x0) - hamiltonian = hamiltonian_function(x0) - eigenvalue, eigenvector = np.linalg.eigh(hamiltonian) - eigenvalue_array[i0, :] = eigenvalue - i0 += 1 - import guan - guan.statistics_of_guan_package() - return eigenvalue_array - -# 输入哈密顿量函数（带两组参数），计算两组参数下的本征值，返回本征值向量组 -def calculate_eigenvalue_with_two_parameters(x_array, y_array, hamiltonian_function, print_show=0, print_show_more=0): - import numpy as np - dim_x = np.array(x_array).shape[0] - dim_y = np.array(y_array).shape[0] - if np.array(hamiltonian_function(0,0)).shape==(): - eigenvalue_array = np.zeros((dim_y, dim_x, 1)) - i0 = 0 - for y0 in y_array: - j0 = 0 - for x0 in x_array: - hamiltonian = hamiltonian_function(x0, y0) - eigenvalue_array[i0, j0, 0] = np.real(hamiltonian) - j0 += 1 - i0 += 1 - else: - dim = np.array(hamiltonian_function(0, 0)).shape[0] - eigenvalue_array = np.zeros((dim_y, dim_x, dim)) - i0 = 0 - for y0 in y_array: - j0 = 0 - if print_show==1: - print(y0) - for x0 in x_array: - if print_show_more==1: - print(x0) - hamiltonian = hamiltonian_function(x0, y0) - eigenvalue, eigenvector = np.linalg.eigh(hamiltonian) - eigenvalue_array[i0, j0, :] = eigenvalue - j0 += 1 - i0 += 1 - import guan - guan.statistics_of_guan_package() - return eigenvalue_array - -# 计算哈密顿量的本征矢 -def calculate_eigenvector(hamiltonian): - import numpy as np - eigenvalue, eigenvector = np.linalg.eigh(hamiltonian) - import guan - guan.statistics_of_guan_package() - return eigenvector - -# 通过二分查找的方法获取和相邻波函数一样规范的波函数 -def find_vector_with_the_same_gauge_with_binary_search(vector_target, vector_ref, show_error=1, show_times=0, show_phase=0, n_test=1000, precision=1e-6): - import numpy as np - import cmath - phase_1_pre = 0 - phase_2_pre = np.pi - for i0 in range(n_test): - test_1 = np.sum(np.abs(vector_target*cmath.exp(1j*phase_1_pre) - vector_ref)) - test_2 = np.sum(np.abs(vector_target*cmath.exp(1j*phase_2_pre) - vector_ref)) - if test_1 < precision: - phase = phase_1_pre - if show_times==1: - print('Binary search times=', i0) - break - if i0 == n_test-1: - phase = phase_1_pre - if show_error==1: - print('Gauge not found with binary search times=', i0) - if test_1 < test_2: - if i0 == 0: - phase_1 = phase_1_pre-(phase_2_pre-phase_1_pre)/2 - phase_2 = phase_1_pre+(phase_2_pre-phase_1_pre)/2 - else: - phase_1 = phase_1_pre - phase_2 = phase_1_pre+(phase_2_pre-phase_1_pre)/2 - else: - if i0 == 0: - phase_1 = phase_2_pre-(phase_2_pre-phase_1_pre)/2 - phase_2 = phase_2_pre+(phase_2_pre-phase_1_pre)/2 - else: - phase_1 = phase_2_pre-(phase_2_pre-phase_1_pre)/2 - phase_2 = phase_2_pre - phase_1_pre = phase_1 - phase_2_pre = phase_2 - vector_target = vector_target*cmath.exp(1j*phase) - if show_phase==1: - print('Phase=', phase) - import guan - guan.statistics_of_guan_package() - return vector_target - -# 通过使得波函数的一个非零分量为实数，得到固定规范的波函数 -def find_vector_with_fixed_gauge_by_making_one_component_real(vector, precision=0.005, index=None): - import numpy as np - import cmath - vector = np.array(vector) - if index == None: - index = np.argmax(np.abs(vector)) - sign_pre = np.sign(np.imag(vector[index])) - for phase in np.arange(0, 2*np.pi, precision): - sign = np.sign(np.imag(vector[index]*cmath.exp(1j*phase))) - if np.abs(np.imag(vector[index]*cmath.exp(1j*phase))) < 1e-9 or sign == -sign_pre: - break - sign_pre = sign - vector = vector*cmath.exp(1j*phase) - if np.real(vector[index]) < 0: - vector = -vector - import guan - guan.statistics_of_guan_package() - return vector - -# 通过使得波函数的一个非零分量为实数，得到固定规范的波函数（在一组波函数中选取最大的那个分量） -def find_vector_array_with_fixed_gauge_by_making_one_component_real(vector_array, precision=0.005): - import numpy as np - import guan - vector_sum = 0 - Num_k = np.array(vector_array).shape[0] - for i0 in range(Num_k): - vector_sum += np.abs(vector_array[i0]) - index = np.argmax(np.abs(vector_sum)) - for i0 in range(Num_k): - vector_array[i0] = guan.find_vector_with_fixed_gauge_by_making_one_component_real(vector_array[i0], precision=precision, index=index) - guan.statistics_of_guan_package() - return vector_array - -# 旋转两个简并的波函数（说明：参数比较多，效率不高） -def rotation_of_degenerate_vectors(vector1, vector2, index1=None, index2=None, precision=0.01, criterion=0.01, show_theta=0): - import numpy as np - import math - import cmath - vector1 = np.array(vector1) - vector2 = np.array(vector2) - if index1 == None: - index1 = np.argmax(np.abs(vector1)) - if index2 == None: - index2 = np.argmax(np.abs(vector2)) - if np.abs(vector1[index2])>criterion or np.abs(vector2[index1])>criterion: - for theta in np.arange(0, 2*math.pi, precision): - if show_theta==1: - print(theta) - for phi1 in np.arange(0, 2*math.pi, precision): - for phi2 in np.arange(0, 2*math.pi, precision): - vector1_test = cmath.exp(1j*phi1)*vector1*math.cos(theta)+cmath.exp(1j*phi2)*vector2*math.sin(theta) - vector2_test = -cmath.exp(-1j*phi2)*vector1*math.sin(theta)+cmath.exp(-1j*phi1)*vector2*math.cos(theta) - if np.abs(vector1_test[index2]) 1.001: - print('Error Alert: scattering matrix is not normalized!') - guan.statistics_of_guan_package() - return transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active - -# 从散射矩阵中，获取散射矩阵的信息 -def information_of_scattering_matrix(transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active): - import numpy as np - if np.array(transmission_matrix).shape==(): - dim = 1 - else: - dim = np.array(transmission_matrix).shape[0] - number_of_active_channels = ind_right_active - number_of_evanescent_channels = dim-ind_right_active - k_of_right_moving_active_channels = np.real(k_right[0:ind_right_active]) - k_of_left_moving_active_channels = np.real(k_left[0:ind_right_active]) - velocity_of_right_moving_active_channels = np.real(velocity_right[0:ind_right_active]) - velocity_of_left_moving_active_channels = np.real(velocity_left[0:ind_right_active]) - transmission_matrix_for_active_channels = np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active])) - reflection_matrix_for_active_channels = np.square(np.abs(reflection_matrix[0:ind_right_active, 0:ind_right_active])) - total_transmission_of_channels = np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active])), axis=0) - total_conductance = np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active]))) - total_reflection_of_channels = np.sum(np.square(np.abs(reflection_matrix[0:ind_right_active, 0:ind_right_active])), axis=0) - sum_of_transmission_and_reflection_of_channels = np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active])), axis=0) + np.sum(np.square(np.abs(reflection_matrix[0:ind_right_active, 0:ind_right_active])), axis=0) - import guan - guan.statistics_of_guan_package() - return number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels - -# 已知h00和h01，计算散射矩阵并获得散射矩阵的信息 -def calculate_scattering_matrix_and_get_information(fermi_energy, h00, h01, length=100): - import guan - transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active = guan.calculate_scattering_matrix(fermi_energy, h00, h01, length=length) - - number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels = guan.information_of_scattering_matrix(transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active) - - guan.statistics_of_guan_package() - - return number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels - -# 从散射矩阵中，打印出散射矩阵的信息 -def print_or_write_scattering_matrix_with_information_of_scattering_matrix(number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels, print_show=1, write_file=0, filename='a', file_format='.txt'): - if print_show == 1: - print('\nActive channel (left or right) = ', number_of_active_channels) - print('Evanescent channel (left or right) = ', number_of_evanescent_channels, '\n') - print('K of right-moving active channels:\n', k_of_right_moving_active_channels) - print('K of left-moving active channels:\n', k_of_left_moving_active_channels, '\n') - print('Velocity of right-moving active channels:\n', velocity_of_right_moving_active_channels) - print('Velocity of left-moving active channels:\n', velocity_of_left_moving_active_channels, '\n') - print('Transmission matrix:\n', transmission_matrix_for_active_channels) - print('Reflection matrix:\n', reflection_matrix_for_active_channels, '\n') - print('Total transmission of channels:\n', total_transmission_of_channels) - print('Total conductance = ', total_conductance, '\n') - print('Total reflection of channels:\n', total_reflection_of_channels) - print('Sum of transmission and reflection of channels:\n', sum_of_transmission_and_reflection_of_channels, '\n') - if write_file == 1: - with open(filename+file_format, 'w') as f: - f.write('Active channel (left or right) = ' + str(number_of_active_channels) + '\n') - f.write('Evanescent channel (left or right) = ' + str(number_of_evanescent_channels) + '\n\n') - f.write('Channel K Velocity\n') - for ind0 in range(number_of_active_channels): - f.write(' '+str(ind0 + 1) + ' | '+str(k_of_right_moving_active_channels[ind0])+' ' + str(velocity_of_right_moving_active_channels[ind0])+'\n') - f.write('\n') - for ind0 in range(number_of_active_channels): - f.write(' -' + str(ind0 + 1) + ' | ' + str(k_of_left_moving_active_channels[ind0]) + ' ' + str(velocity_of_left_moving_active_channels[ind0]) + '\n') - f.write('\nScattering matrix:\n ') - for ind0 in range(number_of_active_channels): - f.write(str(ind0+1)+' ') - f.write('\n') - for ind1 in range(number_of_active_channels): - f.write(' '+str(ind1+1)+' ') - for ind2 in range(number_of_active_channels): - f.write('%f' % transmission_matrix_for_active_channels[ind1, ind2]+' ') - f.write('\n') - f.write('\n') - for ind1 in range(number_of_active_channels): - f.write(' -'+str(ind1+1)+' ') - for ind2 in range(number_of_active_channels): - f.write('%f' % reflection_matrix_for_active_channels[ind1, ind2]+' ') - f.write('\n') - f.write('\n') - f.write('Total transmission of channels:\n'+str(total_transmission_of_channels)+'\n') - f.write('Total conductance = '+str(total_conductance)+'\n') - import guan - guan.statistics_of_guan_package() - -# 已知h00和h01，计算散射矩阵并打印出散射矩阵的信息 -def print_or_write_scattering_matrix(fermi_energy, h00, h01, length=100, print_show=1, write_file=0, filename='a', file_format='.txt'): - import guan - transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active = guan.calculate_scattering_matrix(fermi_energy, h00, h01, length=length) - - number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels = guan.information_of_scattering_matrix(transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active) - - guan.print_or_write_scattering_matrix_with_information_of_scattering_matrix(number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels, print_show=print_show, write_file=write_file, filename=filename, file_format=file_format) - - guan.statistics_of_guan_package() - -# 在无序下，计算散射矩阵 -def calculate_scattering_matrix_with_disorder(fermi_energy, h00, h01, length=100, disorder_intensity=2.0, disorder_concentration=1.0): - import numpy as np - import math - import copy - import guan - h01 = np.array(h01) - if np.array(h00).shape==(): - dim = 1 - else: - dim = np.array(h00).shape[0] - k_right, k_left, velocity_right, velocity_left, f_right, f_left, u_right, u_left, ind_right_active = guan.get_classified_k_velocity_u_and_f(fermi_energy, h00, h01) - right_self_energy = np.dot(h01, f_right) - left_self_energy = np.dot(h01.transpose().conj(), np.linalg.inv(f_left)) - for i0 in range(length): - disorder = np.zeros((dim, dim)) - for dim0 in range(dim): - if np.random.uniform(0, 1)<=disorder_concentration: - disorder[dim0, dim0] = np.random.uniform(-disorder_intensity, disorder_intensity) - if i0 == 0: - green_nn_n = guan.green_function(fermi_energy, h00, broadening=0, self_energy=left_self_energy) - green_00_n = copy.deepcopy(green_nn_n) - green_0n_n = copy.deepcopy(green_nn_n) - green_n0_n = copy.deepcopy(green_nn_n) - elif i0 != length-1: - green_nn_n = guan.green_function_nn_n(fermi_energy, h00+disorder, h01, green_nn_n, broadening=0) - else: - green_nn_n = guan.green_function_nn_n(fermi_energy, h00, h01, green_nn_n, broadening=0, self_energy=right_self_energy) - green_00_n = guan.green_function_ii_n(green_00_n, green_0n_n, h01, green_nn_n, green_n0_n) - green_0n_n = guan.green_function_in_n(green_0n_n, h01, green_nn_n) - green_n0_n = guan.green_function_ni_n(green_nn_n, h01, green_n0_n) - temp = np.dot(h01.transpose().conj(), np.linalg.inv(f_right)-np.linalg.inv(f_left)) - transmission_matrix = np.dot(np.dot(np.linalg.inv(u_right), np.dot(green_n0_n, temp)), u_right) - reflection_matrix = np.dot(np.dot(np.linalg.inv(u_left), np.dot(green_00_n, temp)-np.identity(dim)), u_right) - for dim0 in range(dim): - for dim1 in range(dim): - if_active = guan.if_active_channel(k_right[dim0])*guan.if_active_channel(k_right[dim1]) - if if_active == 1: - transmission_matrix[dim0, dim1] = math.sqrt(np.abs(velocity_right[dim0]/velocity_right[dim1])) * transmission_matrix[dim0, dim1] - reflection_matrix[dim0, dim1] = math.sqrt(np.abs(velocity_left[dim0]/velocity_right[dim1]))*reflection_matrix[dim0, dim1] - else: - transmission_matrix[dim0, dim1] = 0 - reflection_matrix[dim0, dim1] = 0 - sum_of_tran_refl_array = np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active])), axis=0)+np.sum(np.square(np.abs(reflection_matrix[0:ind_right_active, 0:ind_right_active])), axis=0) - for sum_of_tran_refl in sum_of_tran_refl_array: - if sum_of_tran_refl > 1.001: - print('Error Alert: scattering matrix is not normalized!') - guan.statistics_of_guan_package() - return transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active - -# 在无序下，计算散射矩阵，并获取散射矩阵多次计算的平均信息 -def calculate_scattering_matrix_with_disorder_and_get_averaged_information(fermi_energy, h00, h01, length=100, disorder_intensity=2.0, disorder_concentration=1.0, calculation_times=1): - import guan - transmission_matrix_for_active_channels_averaged = 0 - reflection_matrix_for_active_channels_averaged = 0 - for i0 in range(calculation_times): - transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active = guan.calculate_scattering_matrix_with_disorder(fermi_energy, h00, h01, length, disorder_intensity, disorder_concentration) - - number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels = guan.information_of_scattering_matrix(transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active) - - transmission_matrix_for_active_channels_averaged += transmission_matrix_for_active_channels - reflection_matrix_for_active_channels_averaged += reflection_matrix_for_active_channels - transmission_matrix_for_active_channels_averaged = transmission_matrix_for_active_channels_averaged/calculation_times - reflection_matrix_for_active_channels_averaged = reflection_matrix_for_active_channels_averaged/calculation_times - guan.statistics_of_guan_package() - return transmission_matrix_for_active_channels_averaged, reflection_matrix_for_active_channels_averaged - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -# Module 9: topological invariant - -# 通过高效法计算方格子的陈数 -def calculate_chern_number_for_square_lattice_with_efficient_method(hamiltonian_function, precision=100, print_show=0): - import numpy as np - import math - import cmath - import guan - if np.array(hamiltonian_function(0, 0)).shape==(): - dim = 1 - else: - dim = np.array(hamiltonian_function(0, 0)).shape[0] - delta = 2*math.pi/precision - chern_number = np.zeros(dim, dtype=complex) - for kx in np.arange(-math.pi, math.pi, delta): - if print_show == 1: - print(kx) - for ky in np.arange(-math.pi, math.pi, delta): - H = hamiltonian_function(kx, ky) - vector = guan.calculate_eigenvector(H) - H_delta_kx = hamiltonian_function(kx+delta, ky) - vector_delta_kx = guan.calculate_eigenvector(H_delta_kx) - H_delta_ky = hamiltonian_function(kx, ky+delta) - vector_delta_ky = guan.calculate_eigenvector(H_delta_ky) - H_delta_kx_ky = hamiltonian_function(kx+delta, ky+delta) - vector_delta_kx_ky = guan.calculate_eigenvector(H_delta_kx_ky) - for i in range(dim): - vector_i = vector[:, i] - vector_delta_kx_i = vector_delta_kx[:, i] - vector_delta_ky_i = vector_delta_ky[:, i] - vector_delta_kx_ky_i = vector_delta_kx_ky[:, i] - Ux = np.dot(np.conj(vector_i), vector_delta_kx_i)/abs(np.dot(np.conj(vector_i), vector_delta_kx_i)) - Uy = np.dot(np.conj(vector_i), vector_delta_ky_i)/abs(np.dot(np.conj(vector_i), vector_delta_ky_i)) - Ux_y = np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i)) - Uy_x = np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i)) - F = cmath.log(Ux*Uy_x*(1/Ux_y)*(1/Uy)) - chern_number[i] = chern_number[i] + F - chern_number = chern_number/(2*math.pi*1j) - guan.statistics_of_guan_package() - return chern_number - -# 通过高效法计算方格子的陈数（可计算简并的情况） -def calculate_chern_number_for_square_lattice_with_efficient_method_for_degenerate_case(hamiltonian_function, index_of_bands=[0, 1], precision=100, print_show=0): - import numpy as np - import math - import cmath - delta = 2*math.pi/precision - chern_number = 0 - for kx in np.arange(-math.pi, math.pi, delta): - if print_show == 1: - print(kx) - for ky in np.arange(-math.pi, math.pi, delta): - H = hamiltonian_function(kx, ky) - eigenvalue, vector = np.linalg.eigh(H) - H_delta_kx = hamiltonian_function(kx+delta, ky) - eigenvalue, vector_delta_kx = np.linalg.eigh(H_delta_kx) - H_delta_ky = hamiltonian_function(kx, ky+delta) - eigenvalue, vector_delta_ky = np.linalg.eigh(H_delta_ky) - H_delta_kx_ky = hamiltonian_function(kx+delta, ky+delta) - eigenvalue, vector_delta_kx_ky = np.linalg.eigh(H_delta_kx_ky) - dim = len(index_of_bands) - det_value = 1 - # first dot product - dot_matrix = np.zeros((dim , dim), dtype=complex) - i0 = 0 - for dim1 in index_of_bands: - j0 = 0 - for dim2 in index_of_bands: - dot_matrix[i0, j0] = np.dot(np.conj(vector[:, dim1]), vector_delta_kx[:, dim2]) - j0 += 1 - i0 += 1 - dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) - det_value = det_value*dot_matrix - # second dot product - dot_matrix = np.zeros((dim , dim), dtype=complex) - i0 = 0 - for dim1 in index_of_bands: - j0 = 0 - for dim2 in index_of_bands: - dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx[:, dim1]), vector_delta_kx_ky[:, dim2]) - j0 += 1 - i0 += 1 - dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) - det_value = det_value*dot_matrix - # third dot product - dot_matrix = np.zeros((dim , dim), dtype=complex) - i0 = 0 - for dim1 in index_of_bands: - j0 = 0 - for dim2 in index_of_bands: - dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx_ky[:, dim1]), vector_delta_ky[:, dim2]) - j0 += 1 - i0 += 1 - dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) - det_value = det_value*dot_matrix - # four dot product - dot_matrix = np.zeros((dim , dim), dtype=complex) - i0 = 0 - for dim1 in index_of_bands: - j0 = 0 - for dim2 in index_of_bands: - dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_ky[:, dim1]), vector[:, dim2]) - j0 += 1 - i0 += 1 - dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) - det_value= det_value*dot_matrix - chern_number += cmath.log(det_value) - chern_number = chern_number/(2*math.pi*1j) - import guan - guan.statistics_of_guan_package() - return chern_number - -# 通过Wilson loop方法计算方格子的陈数 -def calculate_chern_number_for_square_lattice_with_wilson_loop(hamiltonian_function, precision_of_plaquettes=20, precision_of_wilson_loop=5, print_show=0): - import numpy as np - import math - delta = 2*math.pi/precision_of_plaquettes - chern_number = 0 - for kx in np.arange(-math.pi, math.pi, delta): - if print_show == 1: - print(kx) - for ky in np.arange(-math.pi, math.pi, delta): - vector_array = [] - # line_1 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx+delta/precision_of_wilson_loop*i0, ky) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - # line_2 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx+delta, ky+delta/precision_of_wilson_loop*i0) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - # line_3 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx+delta-delta/precision_of_wilson_loop*i0, ky+delta) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - # line_4 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx, ky+delta-delta/precision_of_wilson_loop*i0) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - wilson_loop = 1 - for i0 in range(len(vector_array)-1): - wilson_loop = wilson_loop*np.dot(vector_array[i0].transpose().conj(), vector_array[i0+1]) - wilson_loop = wilson_loop*np.dot(vector_array[len(vector_array)-1].transpose().conj(), vector_array[0]) - arg = np.log(np.diagonal(wilson_loop))/1j - chern_number = chern_number + arg - chern_number = chern_number/(2*math.pi) - import guan - guan.statistics_of_guan_package() - return chern_number - -# 通过Wilson loop方法计算方格子的陈数（可计算简并的情况） -def calculate_chern_number_for_square_lattice_with_wilson_loop_for_degenerate_case(hamiltonian_function, index_of_bands=[0, 1], precision_of_plaquettes=20, precision_of_wilson_loop=5, print_show=0): - import numpy as np - import math - delta = 2*math.pi/precision_of_plaquettes - chern_number = 0 - for kx in np.arange(-math.pi, math.pi, delta): - if print_show == 1: - print(kx) - for ky in np.arange(-math.pi, math.pi, delta): - vector_array = [] - # line_1 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx+delta/precision_of_wilson_loop*i0, ky) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - # line_2 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx+delta, ky+delta/precision_of_wilson_loop*i0) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - # line_3 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx+delta-delta/precision_of_wilson_loop*i0, ky+delta) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - # line_4 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx, ky+delta-delta/precision_of_wilson_loop*i0) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - wilson_loop = 1 - dim = len(index_of_bands) - for i0 in range(len(vector_array)-1): - dot_matrix = np.zeros((dim , dim), dtype=complex) - i01 = 0 - for dim1 in index_of_bands: - i02 = 0 - for dim2 in index_of_bands: - dot_matrix[i01, i02] = np.dot(vector_array[i0][:, dim1].transpose().conj(), vector_array[i0+1][:, dim2]) - i02 += 1 - i01 += 1 - det_value = np.linalg.det(dot_matrix) - wilson_loop = wilson_loop*det_value - dot_matrix_plus = np.zeros((dim , dim), dtype=complex) - i01 = 0 - for dim1 in index_of_bands: - i02 = 0 - for dim2 in index_of_bands: - dot_matrix_plus[i01, i02] = np.dot(vector_array[len(vector_array)-1][:, dim1].transpose().conj(), vector_array[0][:, dim2]) - i02 += 1 - i01 += 1 - det_value = np.linalg.det(dot_matrix_plus) - wilson_loop = wilson_loop*det_value - arg = np.log(wilson_loop)/1j - chern_number = chern_number + arg - chern_number = chern_number/(2*math.pi) - import guan - guan.statistics_of_guan_package() - return chern_number - -# 通过高效法计算贝利曲率 -def calculate_berry_curvature_with_efficient_method(hamiltonian_function, k_min='default', k_max='default', precision=100, print_show=0): - import numpy as np - import cmath - import guan - import math - if k_min == 'default': - k_min = -math.pi - if k_max == 'default': - k_max=math.pi - if np.array(hamiltonian_function(0, 0)).shape==(): - dim = 1 - else: - dim = np.array(hamiltonian_function(0, 0)).shape[0] - delta = (k_max-k_min)/precision - k_array = np.arange(k_min, k_max, delta) - berry_curvature_array = np.zeros((k_array.shape[0], k_array.shape[0], dim), dtype=complex) - i0 = 0 - for kx in k_array: - if print_show == 1: - print(kx) - j0 = 0 - for ky in k_array: - H = hamiltonian_function(kx, ky) - vector = guan.calculate_eigenvector(H) - H_delta_kx = hamiltonian_function(kx+delta, ky) - vector_delta_kx = guan.calculate_eigenvector(H_delta_kx) - H_delta_ky = hamiltonian_function(kx, ky+delta) - vector_delta_ky = guan.calculate_eigenvector(H_delta_ky) - H_delta_kx_ky = hamiltonian_function(kx+delta, ky+delta) - vector_delta_kx_ky = guan.calculate_eigenvector(H_delta_kx_ky) - for i in range(dim): - vector_i = vector[:, i] - vector_delta_kx_i = vector_delta_kx[:, i] - vector_delta_ky_i = vector_delta_ky[:, i] - vector_delta_kx_ky_i = vector_delta_kx_ky[:, i] - Ux = np.dot(np.conj(vector_i), vector_delta_kx_i)/abs(np.dot(np.conj(vector_i), vector_delta_kx_i)) - Uy = np.dot(np.conj(vector_i), vector_delta_ky_i)/abs(np.dot(np.conj(vector_i), vector_delta_ky_i)) - Ux_y = np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i)) - Uy_x = np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i)) - berry_curvature = cmath.log(Ux*Uy_x*(1/Ux_y)*(1/Uy))/delta/delta*1j - berry_curvature_array[j0, i0, i] = berry_curvature - j0 += 1 - i0 += 1 - guan.statistics_of_guan_package() - return k_array, berry_curvature_array - -# 通过高效法计算贝利曲率（可计算简并的情况） -def calculate_berry_curvature_with_efficient_method_for_degenerate_case(hamiltonian_function, index_of_bands=[0, 1], k_min='default', k_max='default', precision=100, print_show=0): - import numpy as np - import cmath - import math - if k_min == 'default': - k_min = -math.pi - if k_max == 'default': - k_max=math.pi - delta = (k_max-k_min)/precision - k_array = np.arange(k_min, k_max, delta) - berry_curvature_array = np.zeros((k_array.shape[0], k_array.shape[0]), dtype=complex) - i00 = 0 - for kx in np.arange(k_min, k_max, delta): - if print_show == 1: - print(kx) - j00 = 0 - for ky in np.arange(k_min, k_max, delta): - H = hamiltonian_function(kx, ky) - eigenvalue, vector = np.linalg.eigh(H) - H_delta_kx = hamiltonian_function(kx+delta, ky) - eigenvalue, vector_delta_kx = np.linalg.eigh(H_delta_kx) - H_delta_ky = hamiltonian_function(kx, ky+delta) - eigenvalue, vector_delta_ky = np.linalg.eigh(H_delta_ky) - H_delta_kx_ky = hamiltonian_function(kx+delta, ky+delta) - eigenvalue, vector_delta_kx_ky = np.linalg.eigh(H_delta_kx_ky) - dim = len(index_of_bands) - det_value = 1 - # first dot product - dot_matrix = np.zeros((dim , dim), dtype=complex) - i0 = 0 - for dim1 in index_of_bands: - j0 = 0 - for dim2 in index_of_bands: - dot_matrix[i0, j0] = np.dot(np.conj(vector[:, dim1]), vector_delta_kx[:, dim2]) - j0 += 1 - i0 += 1 - dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) - det_value = det_value*dot_matrix - # second dot product - dot_matrix = np.zeros((dim , dim), dtype=complex) - i0 = 0 - for dim1 in index_of_bands: - j0 = 0 - for dim2 in index_of_bands: - dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx[:, dim1]), vector_delta_kx_ky[:, dim2]) - j0 += 1 - i0 += 1 - dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) - det_value = det_value*dot_matrix - # third dot product - dot_matrix = np.zeros((dim , dim), dtype=complex) - i0 = 0 - for dim1 in index_of_bands: - j0 = 0 - for dim2 in index_of_bands: - dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx_ky[:, dim1]), vector_delta_ky[:, dim2]) - j0 += 1 - i0 += 1 - dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) - det_value = det_value*dot_matrix - # four dot product - dot_matrix = np.zeros((dim , dim), dtype=complex) - i0 = 0 - for dim1 in index_of_bands: - j0 = 0 - for dim2 in index_of_bands: - dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_ky[:, dim1]), vector[:, dim2]) - j0 += 1 - i0 += 1 - dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) - det_value= det_value*dot_matrix - berry_curvature = cmath.log(det_value)/delta/delta*1j - berry_curvature_array[j00, i00] = berry_curvature - j00 += 1 - i00 += 1 - import guan - guan.statistics_of_guan_package() - return k_array, berry_curvature_array - -# 通过Wilson loop方法计算贝里曲率 -def calculate_berry_curvature_with_wilson_loop(hamiltonian_function, k_min='default', k_max='default', precision_of_plaquettes=20, precision_of_wilson_loop=5, print_show=0): - import numpy as np - import math - if k_min == 'default': - k_min = -math.pi - if k_max == 'default': - k_max=math.pi - if np.array(hamiltonian_function(0, 0)).shape==(): - dim = 1 - else: - dim = np.array(hamiltonian_function(0, 0)).shape[0] - delta = (k_max-k_min)/precision_of_plaquettes - k_array = np.arange(k_min, k_max, delta) - berry_curvature_array = np.zeros((k_array.shape[0], k_array.shape[0], dim), dtype=complex) - i00 = 0 - for kx in k_array: - if print_show == 1: - print(kx) - j00 = 0 - for ky in k_array: - vector_array = [] - # line_1 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx+delta/precision_of_wilson_loop*i0, ky) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - # line_2 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx+delta, ky+delta/precision_of_wilson_loop*i0) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - # line_3 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx+delta-delta/precision_of_wilson_loop*i0, ky+delta) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - # line_4 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx, ky+delta-delta/precision_of_wilson_loop*i0) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - wilson_loop = 1 - for i0 in range(len(vector_array)-1): - wilson_loop = wilson_loop*np.dot(vector_array[i0].transpose().conj(), vector_array[i0+1]) - wilson_loop = wilson_loop*np.dot(vector_array[len(vector_array)-1].transpose().conj(), vector_array[0]) - berry_curvature = np.log(np.diagonal(wilson_loop))/delta/delta*1j - berry_curvature_array[j00, i00, :]=berry_curvature - j00 += 1 - i00 += 1 - import guan - guan.statistics_of_guan_package() - return k_array, berry_curvature_array - -# 通过Wilson loop方法计算贝里曲率（可计算简并的情况） -def calculate_berry_curvature_with_wilson_loop_for_degenerate_case(hamiltonian_function, index_of_bands=[0, 1], k_min='default', k_max='default', precision_of_plaquettes=20, precision_of_wilson_loop=5, print_show=0): - import numpy as np - import math - if k_min == 'default': - k_min = -math.pi - if k_max == 'default': - k_max=math.pi - delta = (k_max-k_min)/precision_of_plaquettes - k_array = np.arange(k_min, k_max, delta) - berry_curvature_array = np.zeros((k_array.shape[0], k_array.shape[0]), dtype=complex) - i000 = 0 - for kx in k_array: - if print_show == 1: - print(kx) - j000 = 0 - for ky in k_array: - vector_array = [] - # line_1 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx+delta/precision_of_wilson_loop*i0, ky) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - # line_2 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx+delta, ky+delta/precision_of_wilson_loop*i0) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - # line_3 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx+delta-delta/precision_of_wilson_loop*i0, ky+delta) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - # line_4 - for i0 in range(precision_of_wilson_loop): - H_delta = hamiltonian_function(kx, ky+delta-delta/precision_of_wilson_loop*i0) - eigenvalue, eigenvector = np.linalg.eig(H_delta) - vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] - vector_array.append(vector_delta) - wilson_loop = 1 - dim = len(index_of_bands) - for i0 in range(len(vector_array)-1): - dot_matrix = np.zeros((dim , dim), dtype=complex) - i01 = 0 - for dim1 in index_of_bands: - i02 = 0 - for dim2 in index_of_bands: - dot_matrix[i01, i02] = np.dot(vector_array[i0][:, dim1].transpose().conj(), vector_array[i0+1][:, dim2]) - i02 += 1 - i01 += 1 - det_value = np.linalg.det(dot_matrix) - wilson_loop = wilson_loop*det_value - dot_matrix_plus = np.zeros((dim , dim), dtype=complex) - i01 = 0 - for dim1 in index_of_bands: - i02 = 0 - for dim2 in index_of_bands: - dot_matrix_plus[i01, i02] = np.dot(vector_array[len(vector_array)-1][:, dim1].transpose().conj(), vector_array[0][:, dim2]) - i02 += 1 - i01 += 1 - det_value = np.linalg.det(dot_matrix_plus) - wilson_loop = wilson_loop*det_value - berry_curvature = np.log(wilson_loop)/delta/delta*1j - berry_curvature_array[j000, i000]=berry_curvature - j000 += 1 - i000 += 1 - import guan - guan.statistics_of_guan_package() - return k_array, berry_curvature_array - -# 计算蜂窝格子的陈数（高效法） -def calculate_chern_number_for_honeycomb_lattice(hamiltonian_function, a=1, precision=300, print_show=0): - import numpy as np - import math - import cmath - import guan - if np.array(hamiltonian_function(0, 0)).shape==(): - dim = 1 - else: - dim = np.array(hamiltonian_function(0, 0)).shape[0] - chern_number = np.zeros(dim, dtype=complex) - L1 = 4*math.sqrt(3)*math.pi/9/a - L2 = 2*math.sqrt(3)*math.pi/9/a - L3 = 2*math.pi/3/a - delta1 = 2*L1/precision - delta3 = 2*L3/precision - for kx in np.arange(-L1, L1, delta1): - if print_show == 1: - print(kx) - for ky in np.arange(-L3, L3, delta3): - if (-L2<=kx<=L2) or (kx>L2 and -(L1-kx)*math.tan(math.pi/3)<=ky<=(L1-kx)*math.tan(math.pi/3)) or (kx<-L2 and -(kx-(-L1))*math.tan(math.pi/3)<=ky<=(kx-(-L1))*math.tan(math.pi/3)): - H = hamiltonian_function(kx, ky) - vector = guan.calculate_eigenvector(H) - H_delta_kx = hamiltonian_function(kx+delta1, ky) - vector_delta_kx = guan.calculate_eigenvector(H_delta_kx) - H_delta_ky = hamiltonian_function(kx, ky+delta3) - vector_delta_ky = guan.calculate_eigenvector(H_delta_ky) - H_delta_kx_ky = hamiltonian_function(kx+delta1, ky+delta3) - vector_delta_kx_ky = guan.calculate_eigenvector(H_delta_kx_ky) - for i in range(dim): - vector_i = vector[:, i] - vector_delta_kx_i = vector_delta_kx[:, i] - vector_delta_ky_i = vector_delta_ky[:, i] - vector_delta_kx_ky_i = vector_delta_kx_ky[:, i] - Ux = np.dot(np.conj(vector_i), vector_delta_kx_i)/abs(np.dot(np.conj(vector_i), vector_delta_kx_i)) - Uy = np.dot(np.conj(vector_i), vector_delta_ky_i)/abs(np.dot(np.conj(vector_i), vector_delta_ky_i)) - Ux_y = np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i)) - Uy_x = np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i)) - F = cmath.log(Ux*Uy_x*(1/Ux_y)*(1/Uy)) - chern_number[i] = chern_number[i] + F - chern_number = chern_number/(2*math.pi*1j) - guan.statistics_of_guan_package() - return chern_number - -# 计算Wilson loop -def calculate_wilson_loop(hamiltonian_function, k_min='default', k_max='default', precision=100, print_show=0): - import numpy as np - import guan - import math - if k_min == 'default': - k_min = -math.pi - if k_max == 'default': - k_max=math.pi - k_array = np.linspace(k_min, k_max, precision) - dim = np.array(hamiltonian_function(0)).shape[0] - wilson_loop_array = np.ones(dim, dtype=complex) - for i in range(dim): - if print_show == 1: - print(i) - eigenvector_array = [] - for k in k_array: - eigenvector = guan.calculate_eigenvector(hamiltonian_function(k)) - if k != k_max: - eigenvector_array.append(eigenvector[:, i]) - else: - eigenvector_array.append(eigenvector_array[0]) - for i0 in range(precision-1): - F = np.dot(eigenvector_array[i0+1].transpose().conj(), eigenvector_array[i0]) - wilson_loop_array[i] = np.dot(F, wilson_loop_array[i]) - guan.statistics_of_guan_package() - return wilson_loop_array - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -# Module 10: plot figures - -# 导入plt, fig, ax -def import_plt_and_start_fig_ax(adjust_bottom=0.2, adjust_left=0.2, labelsize=20): - import matplotlib.pyplot as plt - fig, ax = plt.subplots() - plt.subplots_adjust(bottom=adjust_bottom, left=adjust_left) - ax.grid() - ax.tick_params(labelsize=labelsize) - labels = ax.get_xticklabels() + ax.get_yticklabels() - [label.set_fontname('Times New Roman') for label in labels] - import guan - guan.statistics_of_guan_package() - return plt, fig, ax - -# 基于plt, fig, ax开始画图 -def plot_without_starting_fig(plt, fig, ax, x_array, y_array, xlabel='x', ylabel='y', title='', fontsize=20, style='', y_min=None, y_max=None, linewidth=None, markersize=None, color=None): - if color==None: - ax.plot(x_array, y_array, style, linewidth=linewidth, markersize=markersize) - else: - ax.plot(x_array, y_array, style, linewidth=linewidth, markersize=markersize, color=color) - ax.set_title(title, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_xlabel(xlabel, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_ylabel(ylabel, fontsize=fontsize, fontfamily='Times New Roman') - if y_min!=None or y_max!=None: - if y_min==None: - y_min=min(y_array) - if y_max==None: - y_max=max(y_array) - ax.set_ylim(y_min, y_max) - import guan - guan.statistics_of_guan_package() - -# 画图 -def plot(x_array, y_array, xlabel='x', ylabel='y', title='', fontsize=20, labelsize=20, show=1, save=0, filename='a', file_format='.jpg', dpi=300, style='', y_min=None, y_max=None, linewidth=None, markersize=None, adjust_bottom=0.2, adjust_left=0.2): - import guan - plt, fig, ax = guan.import_plt_and_start_fig_ax(adjust_bottom=adjust_bottom, adjust_left=adjust_left, labelsize=labelsize) - ax.plot(x_array, y_array, style, linewidth=linewidth, markersize=markersize) - ax.set_title(title, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_xlabel(xlabel, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_ylabel(ylabel, fontsize=fontsize, fontfamily='Times New Roman') - if y_min!=None or y_max!=None: - if y_min==None: - y_min=min(y_array) - if y_max==None: - y_max=max(y_array) - ax.set_ylim(y_min, y_max) - if save == 1: - plt.savefig(filename+file_format, dpi=dpi) - if show == 1: - plt.show() - plt.close('all') - guan.statistics_of_guan_package() - -# 一组横坐标数据，两组纵坐标数据画图 -def plot_two_array(x_array, y1_array, y2_array, xlabel='x', ylabel='y', title='', fontsize=20, labelsize=20, show=1, save=0, filename='a', file_format='.jpg', dpi=300, style_1='', style_2='', y_min=None, y_max=None, linewidth_1=None, linewidth_2=None, markersize_1=None, markersize_2=None, adjust_bottom=0.2, adjust_left=0.2): - import guan - plt, fig, ax = guan.import_plt_and_start_fig_ax(adjust_bottom=adjust_bottom, adjust_left=adjust_left, labelsize=labelsize) - ax.plot(x_array, y1_array, style_1, linewidth=linewidth_1, markersize=markersize_1) - ax.plot(x_array, y2_array, style_2, linewidth=linewidth_2, markersize=markersize_2) - ax.set_title(title, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_xlabel(xlabel, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_ylabel(ylabel, fontsize=fontsize, fontfamily='Times New Roman') - if y_min!=None or y_max!=None: - if y_min==None: - y1_min=min(y1_array) - y2_min=min(y2_array) - y_min=min([y1_min, y2_min]) - if y_max==None: - y1_max=max(y1_array) - y2_max=max(y2_array) - y_max=max([y1_max, y2_max]) - ax.set_ylim(y_min, y_max) - if save == 1: - plt.savefig(filename+file_format, dpi=dpi) - if show == 1: - plt.show() - plt.close('all') - guan.statistics_of_guan_package() - -# 两组横坐标数据，两组纵坐标数据画图 -def plot_two_array_with_two_horizontal_array(x1_array, x2_array, y1_array, y2_array, xlabel='x', ylabel='y', title='', fontsize=20, labelsize=20, show=1, save=0, filename='a', file_format='.jpg', dpi=300, style_1='', style_2='', y_min=None, y_max=None, linewidth_1=None, linewidth_2=None, markersize_1=None, markersize_2=None, adjust_bottom=0.2, adjust_left=0.2): - import guan - plt, fig, ax = guan.import_plt_and_start_fig_ax(adjust_bottom=adjust_bottom, adjust_left=adjust_left, labelsize=labelsize) - ax.plot(x1_array, y1_array, style_1, linewidth=linewidth_1, markersize=markersize_1) - ax.plot(x2_array, y2_array, style_2, linewidth=linewidth_2, markersize=markersize_2) - ax.set_title(title, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_xlabel(xlabel, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_ylabel(ylabel, fontsize=fontsize, fontfamily='Times New Roman') - if y_min!=None or y_max!=None: - if y_min==None: - y1_min=min(y1_array) - y2_min=min(y2_array) - y_min=min([y1_min, y2_min]) - if y_max==None: - y1_max=max(y1_array) - y2_max=max(y2_array) - y_max=max([y1_max, y2_max]) - ax.set_ylim(y_min, y_max) - if save == 1: - plt.savefig(filename+file_format, dpi=dpi) - if show == 1: - plt.show() - plt.close('all') - guan.statistics_of_guan_package() - -# 一组横坐标数据，三组纵坐标数据画图 -def plot_three_array(x_array, y1_array, y2_array, y3_array, xlabel='x', ylabel='y', title='', fontsize=20, labelsize=20, show=1, save=0, filename='a', file_format='.jpg', dpi=300, style_1='', style_2='', style_3='', y_min=None, y_max=None, linewidth_1=None, linewidth_2=None, linewidth_3=None,markersize_1=None, markersize_2=None, markersize_3=None, adjust_bottom=0.2, adjust_left=0.2): - import guan - plt, fig, ax = guan.import_plt_and_start_fig_ax(adjust_bottom=adjust_bottom, adjust_left=adjust_left, labelsize=labelsize) - ax.plot(x_array, y1_array, style_1, linewidth=linewidth_1, markersize=markersize_1) - ax.plot(x_array, y2_array, style_2, linewidth=linewidth_2, markersize=markersize_2) - ax.plot(x_array, y3_array, style_3, linewidth=linewidth_3, markersize=markersize_3) - ax.set_title(title, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_xlabel(xlabel, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_ylabel(ylabel, fontsize=fontsize, fontfamily='Times New Roman') - if y_min!=None or y_max!=None: - if y_min==None: - y1_min=min(y1_array) - y2_min=min(y2_array) - y3_min=min(y3_array) - y_min=min([y1_min, y2_min, y3_min]) - if y_max==None: - y1_max=max(y1_array) - y2_max=max(y2_array) - y3_max=max(y3_array) - y_max=max([y1_max, y2_max, y3_max]) - ax.set_ylim(y_min, y_max) - if save == 1: - plt.savefig(filename+file_format, dpi=dpi) - if show == 1: - plt.show() - plt.close('all') - guan.statistics_of_guan_package() - -# 三组横坐标数据，三组纵坐标数据画图 -def plot_three_array_with_three_horizontal_array(x1_array, x2_array, x3_array, y1_array, y2_array, y3_array, xlabel='x', ylabel='y', title='', fontsize=20, labelsize=20, show=1, save=0, filename='a', file_format='.jpg', dpi=300, style_1='', style_2='', style_3='', y_min=None, y_max=None, linewidth_1=None, linewidth_2=None, linewidth_3=None,markersize_1=None, markersize_2=None, markersize_3=None, adjust_bottom=0.2, adjust_left=0.2): - import guan - plt, fig, ax = guan.import_plt_and_start_fig_ax(adjust_bottom=adjust_bottom, adjust_left=adjust_left, labelsize=labelsize) - ax.plot(x1_array, y1_array, style_1, linewidth=linewidth_1, markersize=markersize_1) - ax.plot(x2_array, y2_array, style_2, linewidth=linewidth_2, markersize=markersize_2) - ax.plot(x3_array, y3_array, style_3, linewidth=linewidth_3, markersize=markersize_3) - ax.set_title(title, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_xlabel(xlabel, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_ylabel(ylabel, fontsize=fontsize, fontfamily='Times New Roman') - if y_min!=None or y_max!=None: - if y_min==None: - y1_min=min(y1_array) - y2_min=min(y2_array) - y3_min=min(y3_array) - y_min=min([y1_min, y2_min, y3_min]) - if y_max==None: - y1_max=max(y1_array) - y2_max=max(y2_array) - y3_max=max(y3_array) - y_max=max([y1_max, y2_max, y3_max]) - ax.set_ylim(y_min, y_max) - if save == 1: - plt.savefig(filename+file_format, dpi=dpi) - if show == 1: - plt.show() - plt.close('all') - guan.statistics_of_guan_package() - -# 画三维图 -def plot_3d_surface(x_array, y_array, matrix, xlabel='x', ylabel='y', zlabel='z', title='', fontsize=20, labelsize=15, show=1, save=0, filename='a', file_format='.jpg', dpi=300, z_min=None, z_max=None, rcount=100, ccount=100): - import numpy as np - import matplotlib.pyplot as plt - from matplotlib import cm - from matplotlib.ticker import LinearLocator - matrix = np.array(matrix) - fig, ax = plt.subplots(subplot_kw={"projection": "3d"}) - plt.subplots_adjust(bottom=0.1, right=0.65) - x_array, y_array = np.meshgrid(x_array, y_array) - if len(matrix.shape) == 2: - surf = ax.plot_surface(x_array, y_array, matrix, rcount=rcount, ccount=ccount, cmap=cm.coolwarm, linewidth=0, antialiased=False) - elif len(matrix.shape) == 3: - for i0 in range(matrix.shape[2]): - surf = ax.plot_surface(x_array, y_array, matrix[:,:,i0], rcount=rcount, ccount=ccount, cmap=cm.coolwarm, linewidth=0, antialiased=False) - ax.set_title(title, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_xlabel(xlabel, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_ylabel(ylabel, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_zlabel(zlabel, fontsize=fontsize, fontfamily='Times New Roman') - ax.zaxis.set_major_locator(LinearLocator(5)) - ax.zaxis.set_major_formatter('{x:.2f}') - if z_min!=None or z_max!=None: - if z_min==None: - z_min=matrix.min() - if z_max==None: - z_max=matrix.max() - ax.set_zlim(z_min, z_max) - ax.tick_params(labelsize=labelsize) - labels = ax.get_xticklabels() + ax.get_yticklabels() + ax.get_zticklabels() - [label.set_fontname('Times New Roman') for label in labels] - cax = plt.axes([0.8, 0.1, 0.05, 0.8]) - cbar = fig.colorbar(surf, cax=cax) - cbar.ax.tick_params(labelsize=labelsize) - for l in cbar.ax.yaxis.get_ticklabels(): - l.set_family('Times New Roman') - if save == 1: - plt.savefig(filename+file_format, dpi=dpi) - if show == 1: - plt.show() - plt.close('all') - import guan - guan.statistics_of_guan_package() - -# 画Contour图 -def plot_contour(x_array, y_array, matrix, xlabel='x', ylabel='y', title='', fontsize=20, labelsize=15, cmap='jet', levels=None, show=1, save=0, filename='a', file_format='.jpg', dpi=300): - import numpy as np - import matplotlib.pyplot as plt - fig, ax = plt.subplots() - plt.subplots_adjust(bottom=0.2, right=0.75, left=0.2) - x_array, y_array = np.meshgrid(x_array, y_array) - contour = ax.contourf(x_array,y_array,matrix,cmap=cmap, levels=levels) - ax.set_title(title, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_xlabel(xlabel, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_ylabel(ylabel, fontsize=fontsize, fontfamily='Times New Roman') - ax.tick_params(labelsize=labelsize) - labels = ax.get_xticklabels() + ax.get_yticklabels() - [label.set_fontname('Times New Roman') for label in labels] - cax = plt.axes([0.8, 0.2, 0.05, 0.68]) - cbar = fig.colorbar(contour, cax=cax) - cbar.ax.tick_params(labelsize=labelsize) - for l in cbar.ax.yaxis.get_ticklabels(): - l.set_family('Times New Roman') - if save == 1: - plt.savefig(filename+file_format, dpi=dpi) - if show == 1: - plt.show() - plt.close('all') - import guan - guan.statistics_of_guan_package() - -# 画棋盘图/伪彩色图 -def plot_pcolor(x_array, y_array, matrix, xlabel='x', ylabel='y', title='', fontsize=20, labelsize=15, cmap='jet', levels=None, show=1, save=0, filename='a', file_format='.jpg', dpi=300): - import numpy as np - import matplotlib.pyplot as plt - fig, ax = plt.subplots() - plt.subplots_adjust(bottom=0.2, right=0.75, left=0.2) - x_array, y_array = np.meshgrid(x_array, y_array) - contour = ax.pcolor(x_array,y_array,matrix, cmap=cmap) - ax.set_title(title, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_xlabel(xlabel, fontsize=fontsize, fontfamily='Times New Roman') - ax.set_ylabel(ylabel, fontsize=fontsize, fontfamily='Times New Roman') - ax.tick_params(labelsize=labelsize) - labels = ax.get_xticklabels() + ax.get_yticklabels() - [label.set_fontname('Times New Roman') for label in labels] - cax = plt.axes([0.8, 0.2, 0.05, 0.68]) - cbar = fig.colorbar(contour, cax=cax) - cbar.ax.tick_params(labelsize=labelsize) - for l in cbar.ax.yaxis.get_ticklabels(): - l.set_family('Times New Roman') - if save == 1: - plt.savefig(filename+file_format, dpi=dpi) - if show == 1: - plt.show() - plt.close('all') - import guan - guan.statistics_of_guan_package() - -# 通过坐标画点和线 -def draw_dots_and_lines(coordinate_array, draw_dots=1, draw_lines=1, max_distance=1.1, line_style='-k', linewidth=1, dot_style='ro', markersize=3, show=1, save=0, filename='a', file_format='.eps', dpi=300): - import numpy as np - import matplotlib.pyplot as plt - coordinate_array = np.array(coordinate_array) - print(coordinate_array.shape) - x_range = max(coordinate_array[:, 0])-min(coordinate_array[:, 0]) - y_range = max(coordinate_array[:, 1])-min(coordinate_array[:, 1]) - fig, ax = plt.subplots(figsize=(6*x_range/y_range,6)) - plt.subplots_adjust(left=0, bottom=0, right=1, top=1) - plt.axis('off') - if draw_lines==1: - for i1 in range(coordinate_array.shape[0]): - for i2 in range(coordinate_array.shape[0]): - if np.sqrt((coordinate_array[i1, 0] - coordinate_array[i2, 0])**2+(coordinate_array[i1, 1] - coordinate_array[i2, 1])**2) < max_distance: - ax.plot([coordinate_array[i1, 0], coordinate_array[i2, 0]], [coordinate_array[i1, 1], coordinate_array[i2, 1]], line_style, linewidth=linewidth) - if draw_dots==1: - for i in range(coordinate_array.shape[0]): - ax.plot(coordinate_array[i, 0], coordinate_array[i, 1], dot_style, markersize=markersize) - if show==1: - plt.show() - if save==1: - if file_format=='.eps': - plt.savefig(filename+file_format) - else: - plt.savefig(filename+file_format, dpi=dpi) - import guan - guan.statistics_of_guan_package() - -# 合并两个图片 -def combine_two_images(image_path_array, figsize=(16,8), show=0, save=1, filename='a', file_format='.jpg', dpi=300): - import numpy as np - num = np.array(image_path_array).shape[0] - if num != 2: - print('Error: The number of images should be two!') - else: - import matplotlib.pyplot as plt - import matplotlib.image as mpimg - fig = plt.figure(figsize=figsize) - plt.subplots_adjust(left=0, right=1, bottom=0, top=1, wspace=0, hspace=0) - ax1 = fig.add_subplot(121) - ax2 = fig.add_subplot(122) - image_1 = mpimg.imread(image_path_array[0]) - image_2 = mpimg.imread(image_path_array[1]) - ax1.imshow(image_1) - ax2.imshow(image_2) - ax1.axis('off') - ax2.axis('off') - if show == 1: - plt.show() - if save == 1: - plt.savefig(filename+file_format, dpi=dpi) - plt.close('all') - import guan - guan.statistics_of_guan_package() - -# 合并三个图片 -def combine_three_images(image_path_array, figsize=(16,5), show=0, save=1, filename='a', file_format='.jpg', dpi=300): - import numpy as np - num = np.array(image_path_array).shape[0] - if num != 3: - print('Error: The number of images should be three!') - else: - import matplotlib.pyplot as plt - import matplotlib.image as mpimg - fig = plt.figure(figsize=figsize) - plt.subplots_adjust(left=0, right=1, bottom=0, top=1, wspace=0, hspace=0) - ax1 = fig.add_subplot(131) - ax2 = fig.add_subplot(132) - ax3 = fig.add_subplot(133) - image_1 = mpimg.imread(image_path_array[0]) - image_2 = mpimg.imread(image_path_array[1]) - image_3 = mpimg.imread(image_path_array[2]) - ax1.imshow(image_1) - ax2.imshow(image_2) - ax3.imshow(image_3) - ax1.axis('off') - ax2.axis('off') - ax3.axis('off') - if show == 1: - plt.show() - if save == 1: - plt.savefig(filename+file_format, dpi=dpi) - plt.close('all') - import guan - guan.statistics_of_guan_package() - -# 合并四个图片 -def combine_four_images(image_path_array, figsize=(16,16), show=0, save=1, filename='a', file_format='.jpg', dpi=300): - import numpy as np - num = np.array(image_path_array).shape[0] - if num != 4: - print('Error: The number of images should be four!') - else: - import matplotlib.pyplot as plt - import matplotlib.image as mpimg - fig = plt.figure(figsize=figsize) - plt.subplots_adjust(left=0, right=1, bottom=0, top=1, wspace=0, hspace=0) - ax1 = fig.add_subplot(221) - ax2 = fig.add_subplot(222) - ax3 = fig.add_subplot(223) - ax4 = fig.add_subplot(224) - image_1 = mpimg.imread(image_path_array[0]) - image_2 = mpimg.imread(image_path_array[1]) - image_3 = mpimg.imread(image_path_array[2]) - image_4 = mpimg.imread(image_path_array[3]) - ax1.imshow(image_1) - ax2.imshow(image_2) - ax3.imshow(image_3) - ax4.imshow(image_4) - ax1.axis('off') - ax2.axis('off') - ax3.axis('off') - ax4.axis('off') - if show == 1: - plt.show() - if save == 1: - plt.savefig(filename+file_format, dpi=dpi) - plt.close('all') - import guan - guan.statistics_of_guan_package() - -# 对于某个目录中的txt文件，批量读取和画图 -def batch_reading_and_plotting(directory, xlabel='x', ylabel='y'): - import re - import os - import guan - for root, dirs, files in os.walk(directory): - for file in files: - if re.search('^txt.', file[::-1]): - filename = file[:-4] - x_array, y_array = guan.read_one_dimensional_data(filename=filename) - guan.plot(x_array, y_array, xlabel=xlabel, ylabel=ylabel, title=filename, show=0, save=1, filename=filename) - guan.statistics_of_guan_package() - -# 制作GIF动画 -def make_gif(image_path_array, filename='a', duration=0.1): - import imageio - images = [] - for image_path in image_path_array: - im = imageio.imread(image_path) - images.append(im) - imageio.mimsave(filename+'.gif', images, 'GIF', duration=duration) - import guan - guan.statistics_of_guan_package() - -# 选取颜色 -def color_matplotlib(): - color_array = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown', 'tab:pink', 'tab:gray', 'tab:olive', 'tab:cyan'] - import guan - guan.statistics_of_guan_package() - return color_array - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -# Module 11: read and write - -# 将数据存到文件 -def dump_data(data, filename, file_format='.txt'): - import pickle - with open(filename+file_format, 'wb') as f: - pickle.dump(data, f) - import guan - guan.statistics_of_guan_package() - -# 从文件中恢复数据到变量 -def load_data(filename, file_format='.txt'): - import pickle - with open(filename+file_format, 'rb') as f: - data = pickle.load(f) - import guan - guan.statistics_of_guan_package() - return data - -# 读取文件中的一维数据（每一行一组x和y） -def read_one_dimensional_data(filename='a', file_format='.txt'): - import numpy as np - f = open(filename+file_format, 'r') - text = f.read() - f.close() - row_list = np.array(text.split('\n')) - dim_column = np.array(row_list[0].split()).shape[0] - x_array = np.array([]) - y_array = np.array([]) - for row in row_list: - column = np.array(row.split()) - if column.shape[0] != 0: - x_array = np.append(x_array, [float(column[0])], axis=0) - y_row = np.zeros(dim_column-1) - for dim0 in range(dim_column-1): - y_row[dim0] = float(column[dim0+1]) - if np.array(y_array).shape[0] == 0: - y_array = [y_row] - else: - y_array = np.append(y_array, [y_row], axis=0) - import guan - guan.statistics_of_guan_package() - return x_array, y_array - -# 读取文件中的一维数据（每一行一组x和y）（支持复数形式） -def read_one_dimensional_complex_data(filename='a', file_format='.txt'): - import numpy as np - f = open(filename+file_format, 'r') - text = f.read() - f.close() - row_list = np.array(text.split('\n')) - dim_column = np.array(row_list[0].split()).shape[0] - x_array = np.array([]) - y_array = np.array([]) - for row in row_list: - column = np.array(row.split()) - if column.shape[0] != 0: - x_array = np.append(x_array, [complex(column[0])], axis=0) - y_row = np.zeros(dim_column-1, dtype=complex) - for dim0 in range(dim_column-1): - y_row[dim0] = complex(column[dim0+1]) - if np.array(y_array).shape[0] == 0: - y_array = [y_row] - else: - y_array = np.append(y_array, [y_row], axis=0) - import guan - guan.statistics_of_guan_package() - return x_array, y_array - -# 读取文件中的二维数据（第一行和列分别为横纵坐标） -def read_two_dimensional_data(filename='a', file_format='.txt'): - import numpy as np - f = open(filename+file_format, 'r') - text = f.read() - f.close() - row_list = np.array(text.split('\n')) - dim_column = np.array(row_list[0].split()).shape[0] - x_array = np.array([]) - y_array = np.array([]) - matrix = np.array([]) - for i0 in range(row_list.shape[0]): - column = np.array(row_list[i0].split()) - if i0 == 0: - x_str = column[1::] - x_array = np.zeros(x_str.shape[0]) - for i00 in range(x_str.shape[0]): - x_array[i00] = float(x_str[i00]) - elif column.shape[0] != 0: - y_array = np.append(y_array, [float(column[0])], axis=0) - matrix_row = np.zeros(dim_column-1) - for dim0 in range(dim_column-1): - matrix_row[dim0] = float(column[dim0+1]) - if np.array(matrix).shape[0] == 0: - matrix = [matrix_row] - else: - matrix = np.append(matrix, [matrix_row], axis=0) - import guan - guan.statistics_of_guan_package() - return x_array, y_array, matrix - -# 读取文件中的二维数据（第一行和列分别为横纵坐标）（支持复数形式） -def read_two_dimensional_complex_data(filename='a', file_format='.txt'): - import numpy as np - f = open(filename+file_format, 'r') - text = f.read() - f.close() - row_list = np.array(text.split('\n')) - dim_column = np.array(row_list[0].split()).shape[0] - x_array = np.array([]) - y_array = np.array([]) - matrix = np.array([]) - for i0 in range(row_list.shape[0]): - column = np.array(row_list[i0].split()) - if i0 == 0: - x_str = column[1::] - x_array = np.zeros(x_str.shape[0], dtype=complex) - for i00 in range(x_str.shape[0]): - x_array[i00] = complex(x_str[i00]) - elif column.shape[0] != 0: - y_array = np.append(y_array, [complex(column[0])], axis=0) - matrix_row = np.zeros(dim_column-1, dtype=complex) - for dim0 in range(dim_column-1): - matrix_row[dim0] = complex(column[dim0+1]) - if np.array(matrix).shape[0] == 0: - matrix = [matrix_row] - else: - matrix = np.append(matrix, [matrix_row], axis=0) - import guan - guan.statistics_of_guan_package() - return x_array, y_array, matrix - -# 读取文件中的二维数据（不包括x和y） -def read_two_dimensional_data_without_xy_array(filename='a', file_format='.txt'): - import numpy as np - matrix = np.loadtxt(filename+file_format) - import guan - guan.statistics_of_guan_package() - return matrix - -# 打开文件用于新增内容 -def open_file(filename='a', file_format='.txt'): - f = open(filename+file_format, 'a', encoding='UTF-8') - import guan - guan.statistics_of_guan_package() - return f - -# 在文件中写入一维数据（每一行一组x和y） -def write_one_dimensional_data(x_array, y_array, filename='a', file_format='.txt'): - import guan - with open(filename+file_format, 'w', encoding='UTF-8') as f: - guan.write_one_dimensional_data_without_opening_file(x_array, y_array, f) - import guan - guan.statistics_of_guan_package() - -# 在文件中写入一维数据（每一行一组x和y）（需要输入文件） -def write_one_dimensional_data_without_opening_file(x_array, y_array, f): - import numpy as np - x_array = np.array(x_array) - y_array = np.array(y_array) - i0 = 0 - for x0 in x_array: - f.write(str(x0)+' ') - if len(y_array.shape) == 1: - f.write(str(y_array[i0])+'\n') - elif len(y_array.shape) == 2: - for j0 in range(y_array.shape[1]): - f.write(str(y_array[i0, j0])+' ') - f.write('\n') - i0 += 1 - import guan - guan.statistics_of_guan_package() - -# 在文件中写入二维数据（第一行和列分别为横纵坐标） -def write_two_dimensional_data(x_array, y_array, matrix, filename='a', file_format='.txt'): - import guan - with open(filename+file_format, 'w', encoding='UTF-8') as f: - guan.write_two_dimensional_data_without_opening_file(x_array, y_array, matrix, f) - guan.statistics_of_guan_package() - -# 在文件中写入二维数据（第一行和列分别为横纵坐标）（需要输入文件） -def write_two_dimensional_data_without_opening_file(x_array, y_array, matrix, f): - import numpy as np - x_array = np.array(x_array) - y_array = np.array(y_array) - matrix = np.array(matrix) - f.write('0 ') - for x0 in x_array: - f.write(str(x0)+' ') - f.write('\n') - i0 = 0 - for y0 in y_array: - f.write(str(y0)) - j0 = 0 - for x0 in x_array: - f.write(' '+str(matrix[i0, j0])+' ') - j0 += 1 - f.write('\n') - i0 += 1 - import guan - guan.statistics_of_guan_package() - -# 在文件中写入二维数据（不包括x和y） -def write_two_dimensional_data_without_xy_array(matrix, filename='a', file_format='.txt'): - import guan - with open(filename+file_format, 'w', encoding='UTF-8') as f: - guan.write_two_dimensional_data_without_xy_array_and_without_opening_file(matrix, f) - guan.statistics_of_guan_package() - -# 在文件中写入二维数据（不包括x和y）（需要输入文件） -def write_two_dimensional_data_without_xy_array_and_without_opening_file(matrix, f): - for row in matrix: - for element in row: - f.write(str(element)+' ') - f.write('\n') - import guan - guan.statistics_of_guan_package() - -# 以显示编号的样式，打印数组 -def print_array_with_index(array, show_index=1, index_type=0): - if show_index==0: - for i0 in array: - print(i0) - else: - if index_type==0: - index = 0 - for i0 in array: - print(index, i0) - index += 1 - else: - index = 0 - for i0 in array: - index += 1 - print(index, i0) - import guan - guan.statistics_of_guan_package() - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -# Module 12: file processing - -# 自动先后运行程序（串行） -def run_programs_sequentially(program_files=['./a.py', './b.py'], execute='python ', show_time=0): - import os - import time - if show_time == 1: - start = time.time() - i0 = 0 - for program_file in program_files: - i0 += 1 - if show_time == 1: - start_0 = time.time() - os.system(execute+program_file) - if show_time == 1: - end_0 = time.time() - print('Running time of program_'+str(i0)+' = '+str((end_0-start_0)/60)+' min') - if show_time == 1: - end = time.time() - print('Total running time = '+str((end-start)/60)+' min') - import guan - guan.statistics_of_guan_package() - -# 如果不存在文件夹，则新建文件夹 -def make_directory(directory='./test'): - import os - if not os.path.exists(directory): - os.makedirs(directory) - import guan - guan.statistics_of_guan_package() - -# 复制一份文件 -def copy_file(file1='./a.txt', file2='./b.txt'): - import shutil - shutil.copy(file1, file2) - import guan - guan.statistics_of_guan_package() - -# 拼接两个PDF文件 -def combine_two_pdf_files(input_file_1='a.pdf', input_file_2='b.pdf', output_file='combined_file.pdf'): - import PyPDF2 - output_pdf = PyPDF2.PdfWriter() - with open(input_file_1, 'rb') as file1: - pdf1 = PyPDF2.PdfReader(file1) - for page in range(len(pdf1.pages)): - output_pdf.add_page(pdf1.pages[page]) - with open(input_file_2, 'rb') as file2: - pdf2 = PyPDF2.PdfReader(file2) - for page in range(len(pdf2.pages)): - output_pdf.add_page(pdf2.pages[page]) - with open(output_file, 'wb') as combined_file: - output_pdf.write(combined_file) - import guan - guan.statistics_of_guan_package() - -# 将PDF文件转成文本 -def pdf_to_text(pdf_path): - from pdfminer.pdfparser import PDFParser, PDFDocument - from pdfminer.pdfinterp import PDFResourceManager, PDFPageInterpreter - from pdfminer.converter import PDFPageAggregator - from pdfminer.layout import LAParams, LTTextBox - from pdfminer.pdfinterp import PDFTextExtractionNotAllowed - import logging - logging.Logger.propagate = False - logging.getLogger().setLevel(logging.ERROR) - praser = PDFParser(open(pdf_path, 'rb')) - doc = PDFDocument() - praser.set_document(doc) - doc.set_parser(praser) - doc.initialize() - if not doc.is_extractable: - raise PDFTextExtractionNotAllowed - else: - rsrcmgr = PDFResourceManager() - laparams = LAParams() - device = PDFPageAggregator(rsrcmgr, laparams=laparams) - interpreter = PDFPageInterpreter(rsrcmgr, device) - content = '' - for page in doc.get_pages(): - interpreter.process_page(page) - layout = device.get_result() - for x in layout: - if isinstance(x, LTTextBox): - content = content + x.get_text().strip() - import guan - guan.statistics_of_guan_package() - return content - -# 获取PDF文件页数 -def get_pdf_page_number(pdf_path): - import PyPDF2 - pdf_file = open(pdf_path, 'rb') - pdf_reader = PyPDF2.PdfReader(pdf_file) - num_pages = len(pdf_reader.pages) - import guan - guan.statistics_of_guan_package() - return num_pages - -# 获取PDF文件指定页面的内容 -def pdf_to_txt_for_a_specific_page(pdf_path, page_num=1): - import PyPDF2 - pdf_file = open(pdf_path, 'rb') - pdf_reader = PyPDF2.PdfReader(pdf_file) - num_pages = len(pdf_reader.pages) - for page_num0 in range(num_pages): - if page_num0 == page_num-1: - page = pdf_reader.pages[page_num0] - page_text = page.extract_text() - pdf_file.close() - import guan - guan.statistics_of_guan_package() - return page_text - -# 获取PDF文献中的链接。例如: link_starting_form='https://doi.org' -def get_links_from_pdf(pdf_path, link_starting_form=''): - import PyPDF2 - import re - pdfReader = PyPDF2.PdfFileReader(pdf_path) - pages = pdfReader.getNumPages() - i0 = 0 - links = [] - for page in range(pages): - pageSliced = pdfReader.getPage(page) - pageObject = pageSliced.getObject() - if '/Annots' in pageObject.keys(): - ann = pageObject['/Annots'] - old = '' - for a in ann: - u = a.getObject() - if '/A' in u.keys(): - if re.search(re.compile('^'+link_starting_form), u['/A']['/URI']): - if u['/A']['/URI'] != old: - links.append(u['/A']['/URI']) - i0 += 1 - old = u['/A']['/URI'] - import guan - guan.statistics_of_guan_package() - return links - -# 通过Sci-Hub网站下载文献 -def download_with_scihub(address=None, num=1): - from bs4 import BeautifulSoup - import re - import requests - import os - if num==1 and address!=None: - address_array = [address] - else: - address_array = [] - for i in range(num): - address = input('\nInput：') - address_array.append(address) - for address in address_array: - r = requests.post('https://sci-hub.st/', data={'request': address}) - print('\nResponse：', r) - print('Address：', r.url) - soup = BeautifulSoup(r.text, features='lxml') - pdf_URL = soup.embed['src'] - # pdf_URL = soup.iframe['src'] # This is a code line of history version which fails to get pdf URL. - if re.search(re.compile('^https:'), pdf_URL): - pass - else: - pdf_URL = 'https:'+pdf_URL - print('PDF address：', pdf_URL) - name = re.search(re.compile('fdp.*?/'),pdf_URL[::-1]).group()[::-1][1::] - print('PDF name：', name) - print('Directory：', os.getcwd()) - print('\nDownloading...') - r = requests.get(pdf_URL, stream=True) - with open(name, 'wb') as f: - for chunk in r.iter_content(chunk_size=32): - f.write(chunk) - print('Completed!\n') - if num != 1: - print('All completed!\n') - import guan - guan.statistics_of_guan_package() - -# 将文件目录结构写入Markdown文件 -def write_file_list_in_markdown(directory='./', filename='a', reverse_positive_or_negative=1, starting_from_h1=None, banned_file_format=[], hide_file_format=None, divided_line=None, show_second_number=None, show_third_number=None): - import os - f = open(filename+'.md', 'w', encoding="utf-8") - filenames1 = os.listdir(directory) - u0 = 0 - for filename1 in filenames1[::reverse_positive_or_negative]: - filename1_with_path = os.path.join(directory,filename1) - if os.path.isfile(filename1_with_path): - if os.path.splitext(filename1)[1] not in banned_file_format: - if hide_file_format == None: - f.write('+ '+str(filename1)+'\n\n') - else: - f.write('+ '+str(os.path.splitext(filename1)[0])+'\n\n') - else: - u0 += 1 - if divided_line != None and u0 != 1: - f.write('--------\n\n') - if starting_from_h1 == None: - f.write('#') - f.write('# '+str(filename1)+'\n\n') - - filenames2 = os.listdir(filename1_with_path) - i0 = 0 - for filename2 in filenames2[::reverse_positive_or_negative]: - filename2_with_path = os.path.join(directory, filename1, filename2) - if os.path.isfile(filename2_with_path): - if os.path.splitext(filename2)[1] not in banned_file_format: - if hide_file_format == None: - f.write('+ '+str(filename2)+'\n\n') - else: - f.write('+ '+str(os.path.splitext(filename2)[0])+'\n\n') - else: - i0 += 1 - if starting_from_h1 == None: - f.write('#') - if show_second_number != None: - f.write('## '+str(i0)+'. '+str(filename2)+'\n\n') - else: - f.write('## '+str(filename2)+'\n\n') - - j0 = 0 - filenames3 = os.listdir(filename2_with_path) - for filename3 in filenames3[::reverse_positive_or_negative]: - filename3_with_path = os.path.join(directory, filename1, filename2, filename3) - if os.path.isfile(filename3_with_path): - if os.path.splitext(filename3)[1] not in banned_file_format: - if hide_file_format == None: - f.write('+ '+str(filename3)+'\n\n') - else: - f.write('+ '+str(os.path.splitext(filename3)[0])+'\n\n') - else: - j0 += 1 - if starting_from_h1 == None: - f.write('#') - if show_third_number != None: - f.write('### ('+str(j0)+') '+str(filename3)+'\n\n') - else: - f.write('### '+str(filename3)+'\n\n') - - filenames4 = os.listdir(filename3_with_path) - for filename4 in filenames4[::reverse_positive_or_negative]: - filename4_with_path = os.path.join(directory, filename1, filename2, filename3, filename4) - if os.path.isfile(filename4_with_path): - if os.path.splitext(filename4)[1] not in banned_file_format: - if hide_file_format == None: - f.write('+ '+str(filename4)+'\n\n') - else: - f.write('+ '+str(os.path.splitext(filename4)[0])+'\n\n') - else: - if starting_from_h1 == None: - f.write('#') - f.write('#### '+str(filename4)+'\n\n') - - filenames5 = os.listdir(filename4_with_path) - for filename5 in filenames5[::reverse_positive_or_negative]: - filename5_with_path = os.path.join(directory, filename1, filename2, filename3, filename4, filename5) - if os.path.isfile(filename5_with_path): - if os.path.splitext(filename5)[1] not in banned_file_format: - if hide_file_format == None: - f.write('+ '+str(filename5)+'\n\n') - else: - f.write('+ '+str(os.path.splitext(filename5)[0])+'\n\n') - else: - if starting_from_h1 == None: - f.write('#') - f.write('##### '+str(filename5)+'\n\n') - - filenames6 = os.listdir(filename5_with_path) - for filename6 in filenames6[::reverse_positive_or_negative]: - filename6_with_path = os.path.join(directory, filename1, filename2, filename3, filename4, filename5, filename6) - if os.path.isfile(filename6_with_path): - if os.path.splitext(filename6)[1] not in banned_file_format: - if hide_file_format == None: - f.write('+ '+str(filename6)+'\n\n') - else: - f.write('+ '+str(os.path.splitext(filename6)[0])+'\n\n') - else: - if starting_from_h1 == None: - f.write('#') - f.write('###### '+str(filename6)+'\n\n') - f.close() - import guan - guan.statistics_of_guan_package() - -# 查找文件名相同的文件 -def find_repeated_file_with_same_filename(directory='./', ignored_directory_with_words=[], ignored_file_with_words=[], num=1000): - import os - from collections import Counter - file_list = [] - for root, dirs, files in os.walk(directory): - for i0 in range(len(files)): - file_list.append(files[i0]) - for word in ignored_directory_with_words: - if word in root: - file_list.remove(files[i0]) - for word in ignored_file_with_words: - if word in files[i0]: - try: - file_list.remove(files[i0]) - except: - pass - count_file = Counter(file_list).most_common(num) - repeated_file = [] - for item in count_file: - if item[1]>1: - repeated_file.append(item) - import guan - guan.statistics_of_guan_package() - return repeated_file - -# 统计各个子文件夹中的文件数量 -def count_file_in_sub_directory(directory='./', smaller_than_num=None): - import os - from collections import Counter - dirs_list = [] - for root, dirs, files in os.walk(directory): - if dirs != []: - for i0 in range(len(dirs)): - dirs_list.append(root+'/'+dirs[i0]) - for sub_dir in dirs_list: - file_list = [] - for root, dirs, files in os.walk(sub_dir): - for i0 in range(len(files)): - file_list.append(files[i0]) - count_file = len(file_list) - if smaller_than_num == None: - print(sub_dir) - print(count_file) - print() - else: - if count_filei0 and abs(a1-a2)

', html, re.S) - if random_on==1: - random.shuffle(contents) - if reverse==1: - contents.reverse() - for content in contents: - soup2 = BeautifulSoup(content, features='lxml') - all_h2 = soup2.find_all('h2') - for h2 in all_h2: - if re.search('\d*. ', h2.get_text()): - word = re.findall('[a-zA-Z].*', h2.get_text(), re.S)[0] - exist = os.path.exists(directory+word+'.mp3') - if not exist: - try: - if re.search(word, html_file): - r = requests.get("https://file.guanjihuan.com/words/"+directory+word+".mp3", stream=True) - with open(directory+word+'.mp3', 'wb') as f: - for chunk in r.iter_content(chunk_size=32): - f.write(chunk) - except: - pass - print(h2.get_text()) - try: - pygame.mixer.init() - track = pygame.mixer.music.load(directory+word+'.mp3') - pygame.mixer.music.play() - if show_link==1: - print('https://www.ldoceonline.com/dictionary/'+word) - except: - pass - translation = re.findall('

.*?

', content, re.S)[0][3:-4] - if show_translation==1: - time.sleep(translation_time) - print(translation) - time.sleep(rest_time) - pygame.mixer.music.stop() - print() - import guan - guan.statistics_of_guan_package() - -# 播放挑选过后的学术单词 -def play_selected_academic_words(reverse=0, random_on=0, bre_or_ame='ame', show_link=1, rest_time=3): - from bs4 import BeautifulSoup - import re - import urllib.request - import requests - import os - import pygame - import time - import ssl - import random - ssl._create_default_https_context = ssl._create_unverified_context - html = urllib.request.urlopen("https://www.guanjihuan.com/archives/24732").read().decode('utf-8') - if bre_or_ame == 'ame': - directory = 'words_mp3_ameProns/' - elif bre_or_ame == 'bre': - directory = 'words_mp3_breProns/' - exist_directory = os.path.exists(directory) - html_file = urllib.request.urlopen("https://file.guanjihuan.com/words/"+directory).read().decode('utf-8') - if exist_directory == 0: - os.makedirs(directory) - soup = BeautifulSoup(html, features='lxml') - contents = re.findall('
  • \d.*?
    • ', html, re.S) - if random_on==1: - random.shuffle(contents) - if reverse==1: - contents.reverse() - for content in contents: - soup2 = BeautifulSoup(content, features='lxml') - all_li = soup2.find_all('li') - for li in all_li: - if re.search('\d*. ', li.get_text()): - word = re.findall('\s[a-zA-Z].*?\s', li.get_text(), re.S)[0][1:-1] - exist = os.path.exists(directory+word+'.mp3') - if not exist: - try: - if re.search(word, html_file): - r = requests.get("https://file.guanjihuan.com/words/"+directory+word+".mp3", stream=True) - with open(directory+word+'.mp3', 'wb') as f: - for chunk in r.iter_content(chunk_size=32): - f.write(chunk) - except: - pass - print(li.get_text()) - try: - pygame.mixer.init() - track = pygame.mixer.music.load(directory+word+'.mp3') - pygame.mixer.music.play() - if show_link==1: - print('https://www.ldoceonline.com/dictionary/'+word) - except: - pass - time.sleep(rest_time) - pygame.mixer.music.stop() - print() - import guan - guan.statistics_of_guan_package() - -# 播放元素周期表上的单词 -def play_element_words(random_on=0, show_translation=1, show_link=1, translation_time=2, rest_time=1): - from bs4 import BeautifulSoup - import re - import urllib.request - import requests - import os - import pygame - import time - import ssl - import random - ssl._create_default_https_context = ssl._create_unverified_context - html = urllib.request.urlopen("https://www.guanjihuan.com/archives/10897").read().decode('utf-8') - directory = 'prons/' - exist_directory = os.path.exists(directory) - html_file = urllib.request.urlopen("https://file.guanjihuan.com/words/periodic_table_of_elements/"+directory).read().decode('utf-8') - if exist_directory == 0: - os.makedirs(directory) - soup = BeautifulSoup(html, features='lxml') - contents = re.findall('

    ', html, re.S) - if random_on==1: - random.shuffle(contents) - for content in contents: - soup2 = BeautifulSoup(content, features='lxml') - all_h2 = soup2.find_all('h2') - for h2 in all_h2: - if re.search('\d*. ', h2.get_text()): - word = re.findall('[a-zA-Z].* \(', h2.get_text(), re.S)[0][:-2] - exist = os.path.exists(directory+word+'.mp3') - if not exist: - try: - if re.search(word, html_file): - r = requests.get("https://file.guanjihuan.com/words/periodic_table_of_elements/prons/"+word+".mp3", stream=True) - with open(directory+word+'.mp3', 'wb') as f: - for chunk in r.iter_content(chunk_size=32): - f.write(chunk) - except: - pass - print(h2.get_text()) - try: - pygame.mixer.init() - track = pygame.mixer.music.load(directory+word+'.mp3') - pygame.mixer.music.play() - if show_link==1: - print('https://www.merriam-webster.com/dictionary/'+word) - except: - pass - translation = re.findall('

    .*?

    ', content, re.S)[0][3:-4] - if show_translation==1: - time.sleep(translation_time) - print(translation) - time.sleep(rest_time) - pygame.mixer.music.stop() - print() - import guan - guan.statistics_of_guan_package() - -# 获取Guan软件包当前模块的所有函数名 -def get_all_function_names_in_current_module(): - import inspect - current_module = inspect.getmodule(inspect.currentframe()) - function_names = [name for name, obj in inspect.getmembers(current_module) if inspect.isfunction(obj)] - import guan - guan.statistics_of_guan_package() - return function_names - -# 统计Guan软件包中的函数数量 -def count_functions_in_current_module(): - import guan - function_names = guan.get_all_function_names_in_current_module() - num_functions = len(function_names) - guan.statistics_of_guan_package() - return num_functions - -# 获取当前函数名 -def get_current_function_name(): - import inspect - current_function_name = inspect.currentframe().f_code.co_name - import guan - guan.statistics_of_guan_package() - return current_function_name - -# 获取调用本函数的函数名 -def get_calling_function_name(layer=1): - import inspect - caller = inspect.stack()[layer] - calling_function_name = caller.function - return calling_function_name - -# 获取当前日期字符串 -def get_date(bar=True): - import datetime - datetime_date = str(datetime.date.today()) - if bar==False: - datetime_date = datetime_date.replace('-', '') - return datetime_date - -# 获取当前时间字符串 -def get_time(): - import datetime - datetime_time = datetime.datetime.now().strftime('%H:%M:%S') - return datetime_time - -# 获取MAC地址 -def get_mac_address(): - import uuid - mac_address = uuid.UUID(int=uuid.getnode()).hex[-12:].upper() - mac_address = '-'.join([mac_address[i:i+2] for i in range(0, 11, 2)]) - return mac_address - -# Guan软件包的使用统计（不涉及到用户的个人数据） -def statistics_of_guan_package(): - try: - import guan - message_calling = guan.get_calling_function_name(layer=3) - if message_calling == '': - import socket - datetime_date = guan.get_date() - datetime_time = guan.get_time() - current_version = guan.get_current_version('guan') - client_socket = socket.socket(socket.AF_INET, socket.SOCK_STREAM) - client_socket.settimeout(0.5) - client_socket.connect(('py.guanjihuan.com', 12345)) - message = guan.get_calling_function_name(layer=2) - send_message = datetime_date + ' ' + datetime_time + ' version_'+current_version+' guan.' + message+'\n' - client_socket.send(send_message.encode()) - client_socket.close() - client_socket2 = socket.socket(socket.AF_INET, socket.SOCK_STREAM) - client_socket2.settimeout(0.5) - client_socket2.connect(('py.guanjihuan.com', 12345)) - mac_address = guan.get_mac_address() - send_mac_address = 'version_'+current_version+' MAC_address: '+mac_address+'\n' - client_socket2.send(send_mac_address.encode()) - client_socket2.close() - except: - pass - -# 获取Python软件包的最新版本 -def get_latest_version(package_name='guan', timeout=2): - import requests - url = f"https://pypi.org/pypi/{package_name}/json" - try: - response = requests.get(url, timeout=timeout) - except: - return None - if response.status_code == 200: - data = response.json() - latest_version = data["info"]["version"] - return latest_version - else: - return None - -# 获取软件包的本机版本 -def get_current_version(package_name='guan'): - import importlib.metadata - try: - current_version = importlib.metadata.version(package_name) - return current_version - except: - return None - -# Guan软件包升级提示 -def notification_of_upgrade(): - try: - import guan - latest_version = guan.get_latest_version(package_name='guan', timeout=2) - current_version = guan.get_current_version('guan') - if latest_version != None and current_version != None: - if latest_version != current_version: - print('提示：您当前使用的版本是 guan-'+current_version+'，目前已经有最新版本 guan-'+latest_version+'。您可以通过以下命令对软件包进行升级：pip install --upgrade guan') - except: - pass - -notification_of_upgrade() \ No newline at end of file +# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about (Ji-Huan Guan, 关济寰). The primary location of this package is on website https://py.guanjihuan.com. The GitHub location of this package is on https://github.com/guanjihuan/py.guanjihuan.com. + +from basic_functions import * +from Fourier_transform import * +from Hamiltonian_of_finite_size_systems import * +from Hamiltonian_of_models_in_reciprocal_space import * +from band_structures_and_wave_functions import * +from Green_functions import * +from density_of_states import * +from quantum_transport import * +from topological_invariant import * +from plot_figures import * +from read_and_write import * +from file_processing import * +from data_processing import * + +import guan +rand_number = guan.get_random_number(start=1, end=10) +if rand_number == 5: + guan.notification_of_upgrade() \ No newline at end of file diff --git a/PyPI/src/guan/band_structures_and_wave_functions.py b/PyPI/src/guan/band_structures_and_wave_functions.py new file mode 100644 index 0000000..7d51264 --- /dev/null +++ b/PyPI/src/guan/band_structures_and_wave_functions.py @@ -0,0 +1,204 @@ +# Module: band_structures_and_wave_functions + +# 计算哈密顿量的本征值 +def calculate_eigenvalue(hamiltonian): + import numpy as np + if np.array(hamiltonian).shape==(): + eigenvalue = np.real(hamiltonian) + else: + eigenvalue, eigenvector = np.linalg.eigh(hamiltonian) + import guan + guan.statistics_of_guan_package() + return eigenvalue + +# 输入哈密顿量函数（带一组参数），计算一组参数下的本征值，返回本征值向量组 +def calculate_eigenvalue_with_one_parameter(x_array, hamiltonian_function, print_show=0): + import numpy as np + dim_x = np.array(x_array).shape[0] + i0 = 0 + if np.array(hamiltonian_function(0)).shape==(): + eigenvalue_array = np.zeros((dim_x, 1)) + for x0 in x_array: + hamiltonian = hamiltonian_function(x0) + eigenvalue_array[i0, 0] = np.real(hamiltonian) + i0 += 1 + else: + dim = np.array(hamiltonian_function(0)).shape[0] + eigenvalue_array = np.zeros((dim_x, dim)) + for x0 in x_array: + if print_show==1: + print(x0) + hamiltonian = hamiltonian_function(x0) + eigenvalue, eigenvector = np.linalg.eigh(hamiltonian) + eigenvalue_array[i0, :] = eigenvalue + i0 += 1 + import guan + guan.statistics_of_guan_package() + return eigenvalue_array + +# 输入哈密顿量函数（带两组参数），计算两组参数下的本征值，返回本征值向量组 +def calculate_eigenvalue_with_two_parameters(x_array, y_array, hamiltonian_function, print_show=0, print_show_more=0): + import numpy as np + dim_x = np.array(x_array).shape[0] + dim_y = np.array(y_array).shape[0] + if np.array(hamiltonian_function(0,0)).shape==(): + eigenvalue_array = np.zeros((dim_y, dim_x, 1)) + i0 = 0 + for y0 in y_array: + j0 = 0 + for x0 in x_array: + hamiltonian = hamiltonian_function(x0, y0) + eigenvalue_array[i0, j0, 0] = np.real(hamiltonian) + j0 += 1 + i0 += 1 + else: + dim = np.array(hamiltonian_function(0, 0)).shape[0] + eigenvalue_array = np.zeros((dim_y, dim_x, dim)) + i0 = 0 + for y0 in y_array: + j0 = 0 + if print_show==1: + print(y0) + for x0 in x_array: + if print_show_more==1: + print(x0) + hamiltonian = hamiltonian_function(x0, y0) + eigenvalue, eigenvector = np.linalg.eigh(hamiltonian) + eigenvalue_array[i0, j0, :] = eigenvalue + j0 += 1 + i0 += 1 + import guan + guan.statistics_of_guan_package() + return eigenvalue_array + +# 计算哈密顿量的本征矢 +def calculate_eigenvector(hamiltonian): + import numpy as np + eigenvalue, eigenvector = np.linalg.eigh(hamiltonian) + import guan + guan.statistics_of_guan_package() + return eigenvector + +# 通过二分查找的方法获取和相邻波函数一样规范的波函数 +def find_vector_with_the_same_gauge_with_binary_search(vector_target, vector_ref, show_error=1, show_times=0, show_phase=0, n_test=1000, precision=1e-6): + import numpy as np + import cmath + phase_1_pre = 0 + phase_2_pre = np.pi + for i0 in range(n_test): + test_1 = np.sum(np.abs(vector_target*cmath.exp(1j*phase_1_pre) - vector_ref)) + test_2 = np.sum(np.abs(vector_target*cmath.exp(1j*phase_2_pre) - vector_ref)) + if test_1 < precision: + phase = phase_1_pre + if show_times==1: + print('Binary search times=', i0) + break + if i0 == n_test-1: + phase = phase_1_pre + if show_error==1: + print('Gauge not found with binary search times=', i0) + if test_1 < test_2: + if i0 == 0: + phase_1 = phase_1_pre-(phase_2_pre-phase_1_pre)/2 + phase_2 = phase_1_pre+(phase_2_pre-phase_1_pre)/2 + else: + phase_1 = phase_1_pre + phase_2 = phase_1_pre+(phase_2_pre-phase_1_pre)/2 + else: + if i0 == 0: + phase_1 = phase_2_pre-(phase_2_pre-phase_1_pre)/2 + phase_2 = phase_2_pre+(phase_2_pre-phase_1_pre)/2 + else: + phase_1 = phase_2_pre-(phase_2_pre-phase_1_pre)/2 + phase_2 = phase_2_pre + phase_1_pre = phase_1 + phase_2_pre = phase_2 + vector_target = vector_target*cmath.exp(1j*phase) + if show_phase==1: + print('Phase=', phase) + import guan + guan.statistics_of_guan_package() + return vector_target + +# 通过使得波函数的一个非零分量为实数，得到固定规范的波函数 +def find_vector_with_fixed_gauge_by_making_one_component_real(vector, precision=0.005, index=None): + import numpy as np + import cmath + vector = np.array(vector) + if index == None: + index = np.argmax(np.abs(vector)) + sign_pre = np.sign(np.imag(vector[index])) + for phase in np.arange(0, 2*np.pi, precision): + sign = np.sign(np.imag(vector[index]*cmath.exp(1j*phase))) + if np.abs(np.imag(vector[index]*cmath.exp(1j*phase))) < 1e-9 or sign == -sign_pre: + break + sign_pre = sign + vector = vector*cmath.exp(1j*phase) + if np.real(vector[index]) < 0: + vector = -vector + import guan + guan.statistics_of_guan_package() + return vector + +# 通过使得波函数的一个非零分量为实数，得到固定规范的波函数（在一组波函数中选取最大的那个分量） +def find_vector_array_with_fixed_gauge_by_making_one_component_real(vector_array, precision=0.005): + import numpy as np + import guan + vector_sum = 0 + Num_k = np.array(vector_array).shape[0] + for i0 in range(Num_k): + vector_sum += np.abs(vector_array[i0]) + index = np.argmax(np.abs(vector_sum)) + for i0 in range(Num_k): + vector_array[i0] = guan.find_vector_with_fixed_gauge_by_making_one_component_real(vector_array[i0], precision=precision, index=index) + guan.statistics_of_guan_package() + return vector_array + +# 旋转两个简并的波函数（说明：参数比较多，效率不高） +def rotation_of_degenerate_vectors(vector1, vector2, index1=None, index2=None, precision=0.01, criterion=0.01, show_theta=0): + import numpy as np + import math + import cmath + vector1 = np.array(vector1) + vector2 = np.array(vector2) + if index1 == None: + index1 = np.argmax(np.abs(vector1)) + if index2 == None: + index2 = np.argmax(np.abs(vector2)) + if np.abs(vector1[index2])>criterion or np.abs(vector2[index1])>criterion: + for theta in np.arange(0, 2*math.pi, precision): + if show_theta==1: + print(theta) + for phi1 in np.arange(0, 2*math.pi, precision): + for phi2 in np.arange(0, 2*math.pi, precision): + vector1_test = cmath.exp(1j*phi1)*vector1*math.cos(theta)+cmath.exp(1j*phi2)*vector2*math.sin(theta) + vector2_test = -cmath.exp(-1j*phi2)*vector1*math.sin(theta)+cmath.exp(-1j*phi1)*vector2*math.cos(theta) + if np.abs(vector1_test[index2])i0 and abs(a1-a2)

    ', html, re.S) + if random_on==1: + random.shuffle(contents) + if reverse==1: + contents.reverse() + for content in contents: + soup2 = BeautifulSoup(content, features='lxml') + all_h2 = soup2.find_all('h2') + for h2 in all_h2: + if re.search('\d*. ', h2.get_text()): + word = re.findall('[a-zA-Z].*', h2.get_text(), re.S)[0] + exist = os.path.exists(directory+word+'.mp3') + if not exist: + try: + if re.search(word, html_file): + r = requests.get("https://file.guanjihuan.com/words/"+directory+word+".mp3", stream=True) + with open(directory+word+'.mp3', 'wb') as f: + for chunk in r.iter_content(chunk_size=32): + f.write(chunk) + except: + pass + print(h2.get_text()) + try: + pygame.mixer.init() + track = pygame.mixer.music.load(directory+word+'.mp3') + pygame.mixer.music.play() + if show_link==1: + print('https://www.ldoceonline.com/dictionary/'+word) + except: + pass + translation = re.findall('

    .*?

    ', content, re.S)[0][3:-4] + if show_translation==1: + time.sleep(translation_time) + print(translation) + time.sleep(rest_time) + pygame.mixer.music.stop() + print() + import guan + guan.statistics_of_guan_package() + +# 播放挑选过后的学术单词 +def play_selected_academic_words(reverse=0, random_on=0, bre_or_ame='ame', show_link=1, rest_time=3): + from bs4 import BeautifulSoup + import re + import urllib.request + import requests + import os + import pygame + import time + import ssl + import random + ssl._create_default_https_context = ssl._create_unverified_context + html = urllib.request.urlopen("https://www.guanjihuan.com/archives/24732").read().decode('utf-8') + if bre_or_ame == 'ame': + directory = 'words_mp3_ameProns/' + elif bre_or_ame == 'bre': + directory = 'words_mp3_breProns/' + exist_directory = os.path.exists(directory) + html_file = urllib.request.urlopen("https://file.guanjihuan.com/words/"+directory).read().decode('utf-8') + if exist_directory == 0: + os.makedirs(directory) + soup = BeautifulSoup(html, features='lxml') + contents = re.findall('
  • \d.*?
    • ', html, re.S) + if random_on==1: + random.shuffle(contents) + if reverse==1: + contents.reverse() + for content in contents: + soup2 = BeautifulSoup(content, features='lxml') + all_li = soup2.find_all('li') + for li in all_li: + if re.search('\d*. ', li.get_text()): + word = re.findall('\s[a-zA-Z].*?\s', li.get_text(), re.S)[0][1:-1] + exist = os.path.exists(directory+word+'.mp3') + if not exist: + try: + if re.search(word, html_file): + r = requests.get("https://file.guanjihuan.com/words/"+directory+word+".mp3", stream=True) + with open(directory+word+'.mp3', 'wb') as f: + for chunk in r.iter_content(chunk_size=32): + f.write(chunk) + except: + pass + print(li.get_text()) + try: + pygame.mixer.init() + track = pygame.mixer.music.load(directory+word+'.mp3') + pygame.mixer.music.play() + if show_link==1: + print('https://www.ldoceonline.com/dictionary/'+word) + except: + pass + time.sleep(rest_time) + pygame.mixer.music.stop() + print() + import guan + guan.statistics_of_guan_package() + +# 播放元素周期表上的单词 +def play_element_words(random_on=0, show_translation=1, show_link=1, translation_time=2, rest_time=1): + from bs4 import BeautifulSoup + import re + import urllib.request + import requests + import os + import pygame + import time + import ssl + import random + ssl._create_default_https_context = ssl._create_unverified_context + html = urllib.request.urlopen("https://www.guanjihuan.com/archives/10897").read().decode('utf-8') + directory = 'prons/' + exist_directory = os.path.exists(directory) + html_file = urllib.request.urlopen("https://file.guanjihuan.com/words/periodic_table_of_elements/"+directory).read().decode('utf-8') + if exist_directory == 0: + os.makedirs(directory) + soup = BeautifulSoup(html, features='lxml') + contents = re.findall('

    ', html, re.S) + if random_on==1: + random.shuffle(contents) + for content in contents: + soup2 = BeautifulSoup(content, features='lxml') + all_h2 = soup2.find_all('h2') + for h2 in all_h2: + if re.search('\d*. ', h2.get_text()): + word = re.findall('[a-zA-Z].* \(', h2.get_text(), re.S)[0][:-2] + exist = os.path.exists(directory+word+'.mp3') + if not exist: + try: + if re.search(word, html_file): + r = requests.get("https://file.guanjihuan.com/words/periodic_table_of_elements/prons/"+word+".mp3", stream=True) + with open(directory+word+'.mp3', 'wb') as f: + for chunk in r.iter_content(chunk_size=32): + f.write(chunk) + except: + pass + print(h2.get_text()) + try: + pygame.mixer.init() + track = pygame.mixer.music.load(directory+word+'.mp3') + pygame.mixer.music.play() + if show_link==1: + print('https://www.merriam-webster.com/dictionary/'+word) + except: + pass + translation = re.findall('

    .*?

    ', content, re.S)[0][3:-4] + if show_translation==1: + time.sleep(translation_time) + print(translation) + time.sleep(rest_time) + pygame.mixer.music.stop() + print() + import guan + guan.statistics_of_guan_package() + +# 获取Guan软件包当前模块的所有函数名 +def get_all_function_names_in_current_module(): + import inspect + current_module = inspect.getmodule(inspect.currentframe()) + function_names = [name for name, obj in inspect.getmembers(current_module) if inspect.isfunction(obj)] + import guan + guan.statistics_of_guan_package() + return function_names + +# 统计Guan软件包中的函数数量 +def count_functions_in_current_module(): + import guan + function_names = guan.get_all_function_names_in_current_module() + num_functions = len(function_names) + guan.statistics_of_guan_package() + return num_functions + +# 获取当前函数名 +def get_current_function_name(): + import inspect + current_function_name = inspect.currentframe().f_code.co_name + import guan + guan.statistics_of_guan_package() + return current_function_name + +# 获取调用本函数的函数名 +def get_calling_function_name(layer=1): + import inspect + caller = inspect.stack()[layer] + calling_function_name = caller.function + return calling_function_name + +# 获取当前日期字符串 +def get_date(bar=True): + import datetime + datetime_date = str(datetime.date.today()) + if bar==False: + datetime_date = datetime_date.replace('-', '') + return datetime_date + +# 获取当前时间字符串 +def get_time(): + import datetime + datetime_time = datetime.datetime.now().strftime('%H:%M:%S') + return datetime_time + +# 获取MAC地址 +def get_mac_address(): + import uuid + mac_address = uuid.UUID(int=uuid.getnode()).hex[-12:].upper() + mac_address = '-'.join([mac_address[i:i+2] for i in range(0, 11, 2)]) + return mac_address + +# Guan软件包的使用统计（不涉及到用户的个人数据） +def statistics_of_guan_package(): + try: + import guan + message_calling = guan.get_calling_function_name(layer=3) + if message_calling == '': + import socket + datetime_date = guan.get_date() + datetime_time = guan.get_time() + current_version = guan.get_current_version('guan') + client_socket = socket.socket(socket.AF_INET, socket.SOCK_STREAM) + client_socket.settimeout(0.5) + client_socket.connect(('py.guanjihuan.com', 12345)) + mac_address = guan.get_mac_address() + message = guan.get_calling_function_name(layer=2) + send_message = datetime_date + ' ' + datetime_time + ' version_'+current_version + ' MAC_address: '+mac_address+' guan.' + message+'\n' + client_socket.send(send_message.encode()) + client_socket.close() + except: + pass + +# 获取Python软件包的最新版本 +def get_latest_version(package_name='guan', timeout=0.5): + import requests + url = f"https://pypi.org/pypi/{package_name}/json" + try: + response = requests.get(url, timeout=timeout) + except: + return None + if response.status_code == 200: + data = response.json() + latest_version = data["info"]["version"] + return latest_version + else: + return None + +# 获取软件包的本机版本 +def get_current_version(package_name='guan'): + import importlib.metadata + try: + current_version = importlib.metadata.version(package_name) + return current_version + except: + return None + +# Guan软件包升级提示 +def notification_of_upgrade(timeout=0.5): + try: + import guan + latest_version = guan.get_latest_version(package_name='guan', timeout=timeout) + current_version = guan.get_current_version('guan') + if latest_version != None and current_version != None: + if latest_version != current_version: + print('提示：您当前使用的版本是 guan-'+current_version+'，目前已经有最新版本 guan-'+latest_version+'。您可以通过以下命令对软件包进行升级：pip install --upgrade guan') + except: + pass \ No newline at end of file diff --git a/PyPI/src/guan/density_of_states.py b/PyPI/src/guan/density_of_states.py new file mode 100644 index 0000000..ed4b48f --- /dev/null +++ b/PyPI/src/guan/density_of_states.py @@ -0,0 +1,164 @@ +# Module: density_of_states + +# 计算体系的总态密度 +def total_density_of_states(fermi_energy, hamiltonian, broadening=0.01): + import numpy as np + import math + import guan + green = guan.green_function(fermi_energy, hamiltonian, broadening) + total_dos = -np.trace(np.imag(green))/math.pi + guan.statistics_of_guan_package() + return total_dos + +# 对于不同费米能，计算体系的总态密度 +def total_density_of_states_with_fermi_energy_array(fermi_energy_array, hamiltonian, broadening=0.01, print_show=0): + import numpy as np + import guan + dim = np.array(fermi_energy_array).shape[0] + total_dos_array = np.zeros(dim) + i0 = 0 + for fermi_energy in fermi_energy_array: + if print_show == 1: + print(fermi_energy) + total_dos_array[i0] = guan.total_density_of_states(fermi_energy, hamiltonian, broadening) + i0 += 1 + guan.statistics_of_guan_package() + return total_dos_array + +# 计算方格子的局域态密度（其中，哈密顿量的维度为：dim_hamiltonian = N1*N2*internal_degree） +def local_density_of_states_for_square_lattice(fermi_energy, hamiltonian, N1, N2, internal_degree=1, broadening=0.01): + import numpy as np + import math + import guan + green = guan.green_function(fermi_energy, hamiltonian, broadening) + local_dos = np.zeros((N2, N1)) + for i1 in range(N1): + for i2 in range(N2): + for i in range(internal_degree): + local_dos[i2, i1] = local_dos[i2, i1]-np.imag(green[i1*N2*internal_degree+i2*internal_degree+i, i1*N2*internal_degree+i2*internal_degree+i])/math.pi + guan.statistics_of_guan_package() + return local_dos + +# 计算立方格子的局域态密度（其中，哈密顿量的维度为：dim_hamiltonian = N1*N2*N3*internal_degree） +def local_density_of_states_for_cubic_lattice(fermi_energy, hamiltonian, N1, N2, N3, internal_degree=1, broadening=0.01): + import numpy as np + import math + import guan + green = guan.green_function(fermi_energy, hamiltonian, broadening) + local_dos = np.zeros((N3, N2, N1)) + for i1 in range(N1): + for i2 in range(N2): + for i3 in range(N3): + for i in range(internal_degree): + local_dos[i3, i2, i1] = local_dos[i3, i2, i1]-np.imag(green[i1*N2*N3*internal_degree+i2*N3*internal_degree+i3*internal_degree+i, i1*N2*N3*internal_degree+i2*N3*internal_degree+i3*internal_degree+i])/math.pi + guan.statistics_of_guan_package() + return local_dos + +# 利用Dyson方程，计算方格子的局域态密度（其中，h00的维度为：dim_h00 = N2*internal_degree） +def local_density_of_states_for_square_lattice_using_dyson_equation(fermi_energy, h00, h01, N2, N1, internal_degree=1, broadening=0.01): + import numpy as np + import math + import guan + local_dos = np.zeros((N2, N1)) + green_11_1 = guan.green_function(fermi_energy, h00, broadening) + for i1 in range(N1): + green_nn_n_minus = green_11_1 + green_in_n_minus = green_11_1 + green_ni_n_minus = green_11_1 + green_ii_n_minus = green_11_1 + for i2_0 in range(i1): + green_nn_n = guan.green_function_nn_n(fermi_energy, h00, h01, green_nn_n_minus, broadening) + green_nn_n_minus = green_nn_n + if i1!=0: + green_in_n_minus = green_nn_n + green_ni_n_minus = green_nn_n + green_ii_n_minus = green_nn_n + for size_0 in range(N1-1-i1): + green_nn_n = guan.green_function_nn_n(fermi_energy, h00, h01, green_nn_n_minus, broadening) + green_nn_n_minus = green_nn_n + green_ii_n = guan.green_function_ii_n(green_ii_n_minus, green_in_n_minus, h01, green_nn_n, green_ni_n_minus) + green_ii_n_minus = green_ii_n + green_in_n = guan.green_function_in_n(green_in_n_minus, h01, green_nn_n) + green_in_n_minus = green_in_n + green_ni_n = guan.green_function_ni_n(green_nn_n, h01, green_ni_n_minus) + green_ni_n_minus = green_ni_n + for i2 in range(N2): + for i in range(internal_degree): + local_dos[i2, i1] = local_dos[i2, i1] - np.imag(green_ii_n_minus[i2*internal_degree+i, i2*internal_degree+i])/math.pi + guan.statistics_of_guan_package() + return local_dos + +# 利用Dyson方程，计算立方格子的局域态密度（其中，h00的维度为：dim_h00 = N2*N3*internal_degree） +def local_density_of_states_for_cubic_lattice_using_dyson_equation(fermi_energy, h00, h01, N3, N2, N1, internal_degree=1, broadening=0.01): + import numpy as np + import math + import guan + local_dos = np.zeros((N3, N2, N1)) + green_11_1 = guan.green_function(fermi_energy, h00, broadening) + for i1 in range(N1): + green_nn_n_minus = green_11_1 + green_in_n_minus = green_11_1 + green_ni_n_minus = green_11_1 + green_ii_n_minus = green_11_1 + for i1_0 in range(i1): + green_nn_n = guan.green_function_nn_n(fermi_energy, h00, h01, green_nn_n_minus, broadening) + green_nn_n_minus = green_nn_n + if i1!=0: + green_in_n_minus = green_nn_n + green_ni_n_minus = green_nn_n + green_ii_n_minus = green_nn_n + for size_0 in range(N1-1-i1): + green_nn_n = guan.green_function_nn_n(fermi_energy, h00, h01, green_nn_n_minus, broadening) + green_nn_n_minus = green_nn_n + green_ii_n = guan.green_function_ii_n(green_ii_n_minus, green_in_n_minus, h01, green_nn_n, green_ni_n_minus) + green_ii_n_minus = green_ii_n + green_in_n = guan.green_function_in_n(green_in_n_minus, h01, green_nn_n) + green_in_n_minus = green_in_n + green_ni_n = guan.green_function_ni_n(green_nn_n, h01, green_ni_n_minus) + green_ni_n_minus = green_ni_n + for i2 in range(N2): + for i3 in range(N3): + for i in range(internal_degree): + local_dos[i3, i2, i1] = local_dos[i3, i2, i1] -np.imag(green_ii_n_minus[i2*N3*internal_degree+i3*internal_degree+i, i2*N3*internal_degree+i3*internal_degree+i])/math.pi + guan.statistics_of_guan_package() + return local_dos + +# 利用Dyson方程，计算方格子条带（考虑了电极自能）的局域态密度（其中，h00的维度为：dim_h00 = N2*internal_degree） +def local_density_of_states_for_square_lattice_with_self_energy_using_dyson_equation(fermi_energy, h00, h01, N2, N1, right_self_energy, left_self_energy, internal_degree=1, broadening=0.01): + import numpy as np + import math + import guan + local_dos = np.zeros((N2, N1)) + green_11_1 = guan.green_function(fermi_energy, h00+left_self_energy, broadening) + for i1 in range(N1): + green_nn_n_minus = green_11_1 + green_in_n_minus = green_11_1 + green_ni_n_minus = green_11_1 + green_ii_n_minus = green_11_1 + for i2_0 in range(i1): + if i2_0 == N1-1-1: + green_nn_n = guan.green_function_nn_n(fermi_energy, h00+right_self_energy, h01, green_nn_n_minus, broadening) + else: + green_nn_n = guan.green_function_nn_n(fermi_energy, h00, h01, green_nn_n_minus, broadening) + green_nn_n_minus = green_nn_n + if i1!=0: + green_in_n_minus = green_nn_n + green_ni_n_minus = green_nn_n + green_ii_n_minus = green_nn_n + for size_0 in range(N1-1-i1): + if size_0 == N1-1-i1-1: + green_nn_n = guan.green_function_nn_n(fermi_energy, h00+right_self_energy, h01, green_nn_n_minus, broadening) + else: + green_nn_n = guan.green_function_nn_n(fermi_energy, h00, h01, green_nn_n_minus, broadening) + green_nn_n_minus = green_nn_n + green_ii_n = guan.green_function_ii_n(green_ii_n_minus, green_in_n_minus, h01, green_nn_n, green_ni_n_minus) + green_ii_n_minus = green_ii_n + green_in_n = guan.green_function_in_n(green_in_n_minus, h01, green_nn_n) + green_in_n_minus = green_in_n + green_ni_n = guan.green_function_ni_n(green_nn_n, h01, green_ni_n_minus) + green_ni_n_minus = green_ni_n + for i2 in range(N2): + for i in range(internal_degree): + local_dos[i2, i1] = local_dos[i2, i1] - np.imag(green_ii_n_minus[i2*internal_degree+i, i2*internal_degree+i])/math.pi + guan.statistics_of_guan_package() + return local_dos diff --git a/PyPI/src/guan/file_processing.py b/PyPI/src/guan/file_processing.py new file mode 100644 index 0000000..6e410d9 --- /dev/null +++ b/PyPI/src/guan/file_processing.py @@ -0,0 +1,501 @@ +# Module: file_processing + +# 自动先后运行程序（串行） +def run_programs_sequentially(program_files=['./a.py', './b.py'], execute='python ', show_time=0): + import os + import time + if show_time == 1: + start = time.time() + i0 = 0 + for program_file in program_files: + i0 += 1 + if show_time == 1: + start_0 = time.time() + os.system(execute+program_file) + if show_time == 1: + end_0 = time.time() + print('Running time of program_'+str(i0)+' = '+str((end_0-start_0)/60)+' min') + if show_time == 1: + end = time.time() + print('Total running time = '+str((end-start)/60)+' min') + import guan + guan.statistics_of_guan_package() + +# 如果不存在文件夹，则新建文件夹 +def make_directory(directory='./test'): + import os + if not os.path.exists(directory): + os.makedirs(directory) + import guan + guan.statistics_of_guan_package() + +# 复制一份文件 +def copy_file(file1='./a.txt', file2='./b.txt'): + import shutil + shutil.copy(file1, file2) + import guan + guan.statistics_of_guan_package() + +# 拼接两个PDF文件 +def combine_two_pdf_files(input_file_1='a.pdf', input_file_2='b.pdf', output_file='combined_file.pdf'): + import PyPDF2 + output_pdf = PyPDF2.PdfWriter() + with open(input_file_1, 'rb') as file1: + pdf1 = PyPDF2.PdfReader(file1) + for page in range(len(pdf1.pages)): + output_pdf.add_page(pdf1.pages[page]) + with open(input_file_2, 'rb') as file2: + pdf2 = PyPDF2.PdfReader(file2) + for page in range(len(pdf2.pages)): + output_pdf.add_page(pdf2.pages[page]) + with open(output_file, 'wb') as combined_file: + output_pdf.write(combined_file) + import guan + guan.statistics_of_guan_package() + +# 将PDF文件转成文本 +def pdf_to_text(pdf_path): + from pdfminer.pdfparser import PDFParser, PDFDocument + from pdfminer.pdfinterp import PDFResourceManager, PDFPageInterpreter + from pdfminer.converter import PDFPageAggregator + from pdfminer.layout import LAParams, LTTextBox + from pdfminer.pdfinterp import PDFTextExtractionNotAllowed + import logging + logging.Logger.propagate = False + logging.getLogger().setLevel(logging.ERROR) + praser = PDFParser(open(pdf_path, 'rb')) + doc = PDFDocument() + praser.set_document(doc) + doc.set_parser(praser) + doc.initialize() + if not doc.is_extractable: + raise PDFTextExtractionNotAllowed + else: + rsrcmgr = PDFResourceManager() + laparams = LAParams() + device = PDFPageAggregator(rsrcmgr, laparams=laparams) + interpreter = PDFPageInterpreter(rsrcmgr, device) + content = '' + for page in doc.get_pages(): + interpreter.process_page(page) + layout = device.get_result() + for x in layout: + if isinstance(x, LTTextBox): + content = content + x.get_text().strip() + import guan + guan.statistics_of_guan_package() + return content + +# 获取PDF文件页数 +def get_pdf_page_number(pdf_path): + import PyPDF2 + pdf_file = open(pdf_path, 'rb') + pdf_reader = PyPDF2.PdfReader(pdf_file) + num_pages = len(pdf_reader.pages) + import guan + guan.statistics_of_guan_package() + return num_pages + +# 获取PDF文件指定页面的内容 +def pdf_to_txt_for_a_specific_page(pdf_path, page_num=1): + import PyPDF2 + pdf_file = open(pdf_path, 'rb') + pdf_reader = PyPDF2.PdfReader(pdf_file) + num_pages = len(pdf_reader.pages) + for page_num0 in range(num_pages): + if page_num0 == page_num-1: + page = pdf_reader.pages[page_num0] + page_text = page.extract_text() + pdf_file.close() + import guan + guan.statistics_of_guan_package() + return page_text + +# 获取PDF文献中的链接。例如: link_starting_form='https://doi.org' +def get_links_from_pdf(pdf_path, link_starting_form=''): + import PyPDF2 + import re + pdfReader = PyPDF2.PdfFileReader(pdf_path) + pages = pdfReader.getNumPages() + i0 = 0 + links = [] + for page in range(pages): + pageSliced = pdfReader.getPage(page) + pageObject = pageSliced.getObject() + if '/Annots' in pageObject.keys(): + ann = pageObject['/Annots'] + old = '' + for a in ann: + u = a.getObject() + if '/A' in u.keys(): + if re.search(re.compile('^'+link_starting_form), u['/A']['/URI']): + if u['/A']['/URI'] != old: + links.append(u['/A']['/URI']) + i0 += 1 + old = u['/A']['/URI'] + import guan + guan.statistics_of_guan_package() + return links + +# 通过Sci-Hub网站下载文献 +def download_with_scihub(address=None, num=1): + from bs4 import BeautifulSoup + import re + import requests + import os + if num==1 and address!=None: + address_array = [address] + else: + address_array = [] + for i in range(num): + address = input('\nInput：') + address_array.append(address) + for address in address_array: + r = requests.post('https://sci-hub.st/', data={'request': address}) + print('\nResponse：', r) + print('Address：', r.url) + soup = BeautifulSoup(r.text, features='lxml') + pdf_URL = soup.embed['src'] + # pdf_URL = soup.iframe['src'] # This is a code line of history version which fails to get pdf URL. + if re.search(re.compile('^https:'), pdf_URL): + pass + else: + pdf_URL = 'https:'+pdf_URL + print('PDF address：', pdf_URL) + name = re.search(re.compile('fdp.*?/'),pdf_URL[::-1]).group()[::-1][1::] + print('PDF name：', name) + print('Directory：', os.getcwd()) + print('\nDownloading...') + r = requests.get(pdf_URL, stream=True) + with open(name, 'wb') as f: + for chunk in r.iter_content(chunk_size=32): + f.write(chunk) + print('Completed!\n') + if num != 1: + print('All completed!\n') + import guan + guan.statistics_of_guan_package() + +# 将文件目录结构写入Markdown文件 +def write_file_list_in_markdown(directory='./', filename='a', reverse_positive_or_negative=1, starting_from_h1=None, banned_file_format=[], hide_file_format=None, divided_line=None, show_second_number=None, show_third_number=None): + import os + f = open(filename+'.md', 'w', encoding="utf-8") + filenames1 = os.listdir(directory) + u0 = 0 + for filename1 in filenames1[::reverse_positive_or_negative]: + filename1_with_path = os.path.join(directory,filename1) + if os.path.isfile(filename1_with_path): + if os.path.splitext(filename1)[1] not in banned_file_format: + if hide_file_format == None: + f.write('+ '+str(filename1)+'\n\n') + else: + f.write('+ '+str(os.path.splitext(filename1)[0])+'\n\n') + else: + u0 += 1 + if divided_line != None and u0 != 1: + f.write('--------\n\n') + if starting_from_h1 == None: + f.write('#') + f.write('# '+str(filename1)+'\n\n') + + filenames2 = os.listdir(filename1_with_path) + i0 = 0 + for filename2 in filenames2[::reverse_positive_or_negative]: + filename2_with_path = os.path.join(directory, filename1, filename2) + if os.path.isfile(filename2_with_path): + if os.path.splitext(filename2)[1] not in banned_file_format: + if hide_file_format == None: + f.write('+ '+str(filename2)+'\n\n') + else: + f.write('+ '+str(os.path.splitext(filename2)[0])+'\n\n') + else: + i0 += 1 + if starting_from_h1 == None: + f.write('#') + if show_second_number != None: + f.write('## '+str(i0)+'. '+str(filename2)+'\n\n') + else: + f.write('## '+str(filename2)+'\n\n') + + j0 = 0 + filenames3 = os.listdir(filename2_with_path) + for filename3 in filenames3[::reverse_positive_or_negative]: + filename3_with_path = os.path.join(directory, filename1, filename2, filename3) + if os.path.isfile(filename3_with_path): + if os.path.splitext(filename3)[1] not in banned_file_format: + if hide_file_format == None: + f.write('+ '+str(filename3)+'\n\n') + else: + f.write('+ '+str(os.path.splitext(filename3)[0])+'\n\n') + else: + j0 += 1 + if starting_from_h1 == None: + f.write('#') + if show_third_number != None: + f.write('### ('+str(j0)+') '+str(filename3)+'\n\n') + else: + f.write('### '+str(filename3)+'\n\n') + + filenames4 = os.listdir(filename3_with_path) + for filename4 in filenames4[::reverse_positive_or_negative]: + filename4_with_path = os.path.join(directory, filename1, filename2, filename3, filename4) + if os.path.isfile(filename4_with_path): + if os.path.splitext(filename4)[1] not in banned_file_format: + if hide_file_format == None: + f.write('+ '+str(filename4)+'\n\n') + else: + f.write('+ '+str(os.path.splitext(filename4)[0])+'\n\n') + else: + if starting_from_h1 == None: + f.write('#') + f.write('#### '+str(filename4)+'\n\n') + + filenames5 = os.listdir(filename4_with_path) + for filename5 in filenames5[::reverse_positive_or_negative]: + filename5_with_path = os.path.join(directory, filename1, filename2, filename3, filename4, filename5) + if os.path.isfile(filename5_with_path): + if os.path.splitext(filename5)[1] not in banned_file_format: + if hide_file_format == None: + f.write('+ '+str(filename5)+'\n\n') + else: + f.write('+ '+str(os.path.splitext(filename5)[0])+'\n\n') + else: + if starting_from_h1 == None: + f.write('#') + f.write('##### '+str(filename5)+'\n\n') + + filenames6 = os.listdir(filename5_with_path) + for filename6 in filenames6[::reverse_positive_or_negative]: + filename6_with_path = os.path.join(directory, filename1, filename2, filename3, filename4, filename5, filename6) + if os.path.isfile(filename6_with_path): + if os.path.splitext(filename6)[1] not in banned_file_format: + if hide_file_format == None: + f.write('+ '+str(filename6)+'\n\n') + else: + f.write('+ '+str(os.path.splitext(filename6)[0])+'\n\n') + else: + if starting_from_h1 == None: + f.write('#') + f.write('###### '+str(filename6)+'\n\n') + f.close() + import guan + guan.statistics_of_guan_package() + +# 查找文件名相同的文件 +def find_repeated_file_with_same_filename(directory='./', ignored_directory_with_words=[], ignored_file_with_words=[], num=1000): + import os + from collections import Counter + file_list = [] + for root, dirs, files in os.walk(directory): + for i0 in range(len(files)): + file_list.append(files[i0]) + for word in ignored_directory_with_words: + if word in root: + file_list.remove(files[i0]) + for word in ignored_file_with_words: + if word in files[i0]: + try: + file_list.remove(files[i0]) + except: + pass + count_file = Counter(file_list).most_common(num) + repeated_file = [] + for item in count_file: + if item[1]>1: + repeated_file.append(item) + import guan + guan.statistics_of_guan_package() + return repeated_file + +# 统计各个子文件夹中的文件数量 +def count_file_in_sub_directory(directory='./', smaller_than_num=None): + import os + from collections import Counter + dirs_list = [] + for root, dirs, files in os.walk(directory): + if dirs != []: + for i0 in range(len(dirs)): + dirs_list.append(root+'/'+dirs[i0]) + for sub_dir in dirs_list: + file_list = [] + for root, dirs, files in os.walk(sub_dir): + for i0 in range(len(files)): + file_list.append(files[i0]) + count_file = len(file_list) + if smaller_than_num == None: + print(sub_dir) + print(count_file) + print() + else: + if count_file 1.001: + print('Error Alert: scattering matrix is not normalized!') + guan.statistics_of_guan_package() + return transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active + +# 从散射矩阵中，获取散射矩阵的信息 +def information_of_scattering_matrix(transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active): + import numpy as np + if np.array(transmission_matrix).shape==(): + dim = 1 + else: + dim = np.array(transmission_matrix).shape[0] + number_of_active_channels = ind_right_active + number_of_evanescent_channels = dim-ind_right_active + k_of_right_moving_active_channels = np.real(k_right[0:ind_right_active]) + k_of_left_moving_active_channels = np.real(k_left[0:ind_right_active]) + velocity_of_right_moving_active_channels = np.real(velocity_right[0:ind_right_active]) + velocity_of_left_moving_active_channels = np.real(velocity_left[0:ind_right_active]) + transmission_matrix_for_active_channels = np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active])) + reflection_matrix_for_active_channels = np.square(np.abs(reflection_matrix[0:ind_right_active, 0:ind_right_active])) + total_transmission_of_channels = np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active])), axis=0) + total_conductance = np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active]))) + total_reflection_of_channels = np.sum(np.square(np.abs(reflection_matrix[0:ind_right_active, 0:ind_right_active])), axis=0) + sum_of_transmission_and_reflection_of_channels = np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active])), axis=0) + np.sum(np.square(np.abs(reflection_matrix[0:ind_right_active, 0:ind_right_active])), axis=0) + import guan + guan.statistics_of_guan_package() + return number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels + +# 已知h00和h01，计算散射矩阵并获得散射矩阵的信息 +def calculate_scattering_matrix_and_get_information(fermi_energy, h00, h01, length=100): + import guan + transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active = guan.calculate_scattering_matrix(fermi_energy, h00, h01, length=length) + + number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels = guan.information_of_scattering_matrix(transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active) + + guan.statistics_of_guan_package() + + return number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels + +# 从散射矩阵中，打印出散射矩阵的信息 +def print_or_write_scattering_matrix_with_information_of_scattering_matrix(number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels, print_show=1, write_file=0, filename='a', file_format='.txt'): + if print_show == 1: + print('\nActive channel (left or right) = ', number_of_active_channels) + print('Evanescent channel (left or right) = ', number_of_evanescent_channels, '\n') + print('K of right-moving active channels:\n', k_of_right_moving_active_channels) + print('K of left-moving active channels:\n', k_of_left_moving_active_channels, '\n') + print('Velocity of right-moving active channels:\n', velocity_of_right_moving_active_channels) + print('Velocity of left-moving active channels:\n', velocity_of_left_moving_active_channels, '\n') + print('Transmission matrix:\n', transmission_matrix_for_active_channels) + print('Reflection matrix:\n', reflection_matrix_for_active_channels, '\n') + print('Total transmission of channels:\n', total_transmission_of_channels) + print('Total conductance = ', total_conductance, '\n') + print('Total reflection of channels:\n', total_reflection_of_channels) + print('Sum of transmission and reflection of channels:\n', sum_of_transmission_and_reflection_of_channels, '\n') + if write_file == 1: + with open(filename+file_format, 'w') as f: + f.write('Active channel (left or right) = ' + str(number_of_active_channels) + '\n') + f.write('Evanescent channel (left or right) = ' + str(number_of_evanescent_channels) + '\n\n') + f.write('Channel K Velocity\n') + for ind0 in range(number_of_active_channels): + f.write(' '+str(ind0 + 1) + ' | '+str(k_of_right_moving_active_channels[ind0])+' ' + str(velocity_of_right_moving_active_channels[ind0])+'\n') + f.write('\n') + for ind0 in range(number_of_active_channels): + f.write(' -' + str(ind0 + 1) + ' | ' + str(k_of_left_moving_active_channels[ind0]) + ' ' + str(velocity_of_left_moving_active_channels[ind0]) + '\n') + f.write('\nScattering matrix:\n ') + for ind0 in range(number_of_active_channels): + f.write(str(ind0+1)+' ') + f.write('\n') + for ind1 in range(number_of_active_channels): + f.write(' '+str(ind1+1)+' ') + for ind2 in range(number_of_active_channels): + f.write('%f' % transmission_matrix_for_active_channels[ind1, ind2]+' ') + f.write('\n') + f.write('\n') + for ind1 in range(number_of_active_channels): + f.write(' -'+str(ind1+1)+' ') + for ind2 in range(number_of_active_channels): + f.write('%f' % reflection_matrix_for_active_channels[ind1, ind2]+' ') + f.write('\n') + f.write('\n') + f.write('Total transmission of channels:\n'+str(total_transmission_of_channels)+'\n') + f.write('Total conductance = '+str(total_conductance)+'\n') + import guan + guan.statistics_of_guan_package() + +# 已知h00和h01，计算散射矩阵并打印出散射矩阵的信息 +def print_or_write_scattering_matrix(fermi_energy, h00, h01, length=100, print_show=1, write_file=0, filename='a', file_format='.txt'): + import guan + transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active = guan.calculate_scattering_matrix(fermi_energy, h00, h01, length=length) + + number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels = guan.information_of_scattering_matrix(transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active) + + guan.print_or_write_scattering_matrix_with_information_of_scattering_matrix(number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels, print_show=print_show, write_file=write_file, filename=filename, file_format=file_format) + + guan.statistics_of_guan_package() + +# 在无序下，计算散射矩阵 +def calculate_scattering_matrix_with_disorder(fermi_energy, h00, h01, length=100, disorder_intensity=2.0, disorder_concentration=1.0): + import numpy as np + import math + import copy + import guan + h01 = np.array(h01) + if np.array(h00).shape==(): + dim = 1 + else: + dim = np.array(h00).shape[0] + k_right, k_left, velocity_right, velocity_left, f_right, f_left, u_right, u_left, ind_right_active = guan.get_classified_k_velocity_u_and_f(fermi_energy, h00, h01) + right_self_energy = np.dot(h01, f_right) + left_self_energy = np.dot(h01.transpose().conj(), np.linalg.inv(f_left)) + for i0 in range(length): + disorder = np.zeros((dim, dim)) + for dim0 in range(dim): + if np.random.uniform(0, 1)<=disorder_concentration: + disorder[dim0, dim0] = np.random.uniform(-disorder_intensity, disorder_intensity) + if i0 == 0: + green_nn_n = guan.green_function(fermi_energy, h00, broadening=0, self_energy=left_self_energy) + green_00_n = copy.deepcopy(green_nn_n) + green_0n_n = copy.deepcopy(green_nn_n) + green_n0_n = copy.deepcopy(green_nn_n) + elif i0 != length-1: + green_nn_n = guan.green_function_nn_n(fermi_energy, h00+disorder, h01, green_nn_n, broadening=0) + else: + green_nn_n = guan.green_function_nn_n(fermi_energy, h00, h01, green_nn_n, broadening=0, self_energy=right_self_energy) + green_00_n = guan.green_function_ii_n(green_00_n, green_0n_n, h01, green_nn_n, green_n0_n) + green_0n_n = guan.green_function_in_n(green_0n_n, h01, green_nn_n) + green_n0_n = guan.green_function_ni_n(green_nn_n, h01, green_n0_n) + temp = np.dot(h01.transpose().conj(), np.linalg.inv(f_right)-np.linalg.inv(f_left)) + transmission_matrix = np.dot(np.dot(np.linalg.inv(u_right), np.dot(green_n0_n, temp)), u_right) + reflection_matrix = np.dot(np.dot(np.linalg.inv(u_left), np.dot(green_00_n, temp)-np.identity(dim)), u_right) + for dim0 in range(dim): + for dim1 in range(dim): + if_active = guan.if_active_channel(k_right[dim0])*guan.if_active_channel(k_right[dim1]) + if if_active == 1: + transmission_matrix[dim0, dim1] = math.sqrt(np.abs(velocity_right[dim0]/velocity_right[dim1])) * transmission_matrix[dim0, dim1] + reflection_matrix[dim0, dim1] = math.sqrt(np.abs(velocity_left[dim0]/velocity_right[dim1]))*reflection_matrix[dim0, dim1] + else: + transmission_matrix[dim0, dim1] = 0 + reflection_matrix[dim0, dim1] = 0 + sum_of_tran_refl_array = np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active])), axis=0)+np.sum(np.square(np.abs(reflection_matrix[0:ind_right_active, 0:ind_right_active])), axis=0) + for sum_of_tran_refl in sum_of_tran_refl_array: + if sum_of_tran_refl > 1.001: + print('Error Alert: scattering matrix is not normalized!') + guan.statistics_of_guan_package() + return transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active + +# 在无序下，计算散射矩阵，并获取散射矩阵多次计算的平均信息 +def calculate_scattering_matrix_with_disorder_and_get_averaged_information(fermi_energy, h00, h01, length=100, disorder_intensity=2.0, disorder_concentration=1.0, calculation_times=1): + import guan + transmission_matrix_for_active_channels_averaged = 0 + reflection_matrix_for_active_channels_averaged = 0 + for i0 in range(calculation_times): + transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active = guan.calculate_scattering_matrix_with_disorder(fermi_energy, h00, h01, length, disorder_intensity, disorder_concentration) + + number_of_active_channels, number_of_evanescent_channels, k_of_right_moving_active_channels, k_of_left_moving_active_channels, velocity_of_right_moving_active_channels, velocity_of_left_moving_active_channels, transmission_matrix_for_active_channels, reflection_matrix_for_active_channels, total_transmission_of_channels, total_conductance, total_reflection_of_channels, sum_of_transmission_and_reflection_of_channels = guan.information_of_scattering_matrix(transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active) + + transmission_matrix_for_active_channels_averaged += transmission_matrix_for_active_channels + reflection_matrix_for_active_channels_averaged += reflection_matrix_for_active_channels + transmission_matrix_for_active_channels_averaged = transmission_matrix_for_active_channels_averaged/calculation_times + reflection_matrix_for_active_channels_averaged = reflection_matrix_for_active_channels_averaged/calculation_times + guan.statistics_of_guan_package() + return transmission_matrix_for_active_channels_averaged, reflection_matrix_for_active_channels_averaged \ No newline at end of file diff --git a/PyPI/src/guan/read_and_write.py b/PyPI/src/guan/read_and_write.py new file mode 100644 index 0000000..cd1e308 --- /dev/null +++ b/PyPI/src/guan/read_and_write.py @@ -0,0 +1,235 @@ +# Module: read_and_write + +# 将数据存到文件 +def dump_data(data, filename, file_format='.txt'): + import pickle + with open(filename+file_format, 'wb') as f: + pickle.dump(data, f) + import guan + guan.statistics_of_guan_package() + +# 从文件中恢复数据到变量 +def load_data(filename, file_format='.txt'): + import pickle + with open(filename+file_format, 'rb') as f: + data = pickle.load(f) + import guan + guan.statistics_of_guan_package() + return data + +# 读取文件中的一维数据（每一行一组x和y） +def read_one_dimensional_data(filename='a', file_format='.txt'): + import numpy as np + f = open(filename+file_format, 'r') + text = f.read() + f.close() + row_list = np.array(text.split('\n')) + dim_column = np.array(row_list[0].split()).shape[0] + x_array = np.array([]) + y_array = np.array([]) + for row in row_list: + column = np.array(row.split()) + if column.shape[0] != 0: + x_array = np.append(x_array, [float(column[0])], axis=0) + y_row = np.zeros(dim_column-1) + for dim0 in range(dim_column-1): + y_row[dim0] = float(column[dim0+1]) + if np.array(y_array).shape[0] == 0: + y_array = [y_row] + else: + y_array = np.append(y_array, [y_row], axis=0) + import guan + guan.statistics_of_guan_package() + return x_array, y_array + +# 读取文件中的一维数据（每一行一组x和y）（支持复数形式） +def read_one_dimensional_complex_data(filename='a', file_format='.txt'): + import numpy as np + f = open(filename+file_format, 'r') + text = f.read() + f.close() + row_list = np.array(text.split('\n')) + dim_column = np.array(row_list[0].split()).shape[0] + x_array = np.array([]) + y_array = np.array([]) + for row in row_list: + column = np.array(row.split()) + if column.shape[0] != 0: + x_array = np.append(x_array, [complex(column[0])], axis=0) + y_row = np.zeros(dim_column-1, dtype=complex) + for dim0 in range(dim_column-1): + y_row[dim0] = complex(column[dim0+1]) + if np.array(y_array).shape[0] == 0: + y_array = [y_row] + else: + y_array = np.append(y_array, [y_row], axis=0) + import guan + guan.statistics_of_guan_package() + return x_array, y_array + +# 读取文件中的二维数据（第一行和列分别为横纵坐标） +def read_two_dimensional_data(filename='a', file_format='.txt'): + import numpy as np + f = open(filename+file_format, 'r') + text = f.read() + f.close() + row_list = np.array(text.split('\n')) + dim_column = np.array(row_list[0].split()).shape[0] + x_array = np.array([]) + y_array = np.array([]) + matrix = np.array([]) + for i0 in range(row_list.shape[0]): + column = np.array(row_list[i0].split()) + if i0 == 0: + x_str = column[1::] + x_array = np.zeros(x_str.shape[0]) + for i00 in range(x_str.shape[0]): + x_array[i00] = float(x_str[i00]) + elif column.shape[0] != 0: + y_array = np.append(y_array, [float(column[0])], axis=0) + matrix_row = np.zeros(dim_column-1) + for dim0 in range(dim_column-1): + matrix_row[dim0] = float(column[dim0+1]) + if np.array(matrix).shape[0] == 0: + matrix = [matrix_row] + else: + matrix = np.append(matrix, [matrix_row], axis=0) + import guan + guan.statistics_of_guan_package() + return x_array, y_array, matrix + +# 读取文件中的二维数据（第一行和列分别为横纵坐标）（支持复数形式） +def read_two_dimensional_complex_data(filename='a', file_format='.txt'): + import numpy as np + f = open(filename+file_format, 'r') + text = f.read() + f.close() + row_list = np.array(text.split('\n')) + dim_column = np.array(row_list[0].split()).shape[0] + x_array = np.array([]) + y_array = np.array([]) + matrix = np.array([]) + for i0 in range(row_list.shape[0]): + column = np.array(row_list[i0].split()) + if i0 == 0: + x_str = column[1::] + x_array = np.zeros(x_str.shape[0], dtype=complex) + for i00 in range(x_str.shape[0]): + x_array[i00] = complex(x_str[i00]) + elif column.shape[0] != 0: + y_array = np.append(y_array, [complex(column[0])], axis=0) + matrix_row = np.zeros(dim_column-1, dtype=complex) + for dim0 in range(dim_column-1): + matrix_row[dim0] = complex(column[dim0+1]) + if np.array(matrix).shape[0] == 0: + matrix = [matrix_row] + else: + matrix = np.append(matrix, [matrix_row], axis=0) + import guan + guan.statistics_of_guan_package() + return x_array, y_array, matrix + +# 读取文件中的二维数据（不包括x和y） +def read_two_dimensional_data_without_xy_array(filename='a', file_format='.txt'): + import numpy as np + matrix = np.loadtxt(filename+file_format) + import guan + guan.statistics_of_guan_package() + return matrix + +# 打开文件用于新增内容 +def open_file(filename='a', file_format='.txt'): + f = open(filename+file_format, 'a', encoding='UTF-8') + import guan + guan.statistics_of_guan_package() + return f + +# 在文件中写入一维数据（每一行一组x和y） +def write_one_dimensional_data(x_array, y_array, filename='a', file_format='.txt'): + import guan + with open(filename+file_format, 'w', encoding='UTF-8') as f: + guan.write_one_dimensional_data_without_opening_file(x_array, y_array, f) + import guan + guan.statistics_of_guan_package() + +# 在文件中写入一维数据（每一行一组x和y）（需要输入文件） +def write_one_dimensional_data_without_opening_file(x_array, y_array, f): + import numpy as np + x_array = np.array(x_array) + y_array = np.array(y_array) + i0 = 0 + for x0 in x_array: + f.write(str(x0)+' ') + if len(y_array.shape) == 1: + f.write(str(y_array[i0])+'\n') + elif len(y_array.shape) == 2: + for j0 in range(y_array.shape[1]): + f.write(str(y_array[i0, j0])+' ') + f.write('\n') + i0 += 1 + import guan + guan.statistics_of_guan_package() + +# 在文件中写入二维数据（第一行和列分别为横纵坐标） +def write_two_dimensional_data(x_array, y_array, matrix, filename='a', file_format='.txt'): + import guan + with open(filename+file_format, 'w', encoding='UTF-8') as f: + guan.write_two_dimensional_data_without_opening_file(x_array, y_array, matrix, f) + guan.statistics_of_guan_package() + +# 在文件中写入二维数据（第一行和列分别为横纵坐标）（需要输入文件） +def write_two_dimensional_data_without_opening_file(x_array, y_array, matrix, f): + import numpy as np + x_array = np.array(x_array) + y_array = np.array(y_array) + matrix = np.array(matrix) + f.write('0 ') + for x0 in x_array: + f.write(str(x0)+' ') + f.write('\n') + i0 = 0 + for y0 in y_array: + f.write(str(y0)) + j0 = 0 + for x0 in x_array: + f.write(' '+str(matrix[i0, j0])+' ') + j0 += 1 + f.write('\n') + i0 += 1 + import guan + guan.statistics_of_guan_package() + +# 在文件中写入二维数据（不包括x和y） +def write_two_dimensional_data_without_xy_array(matrix, filename='a', file_format='.txt'): + import guan + with open(filename+file_format, 'w', encoding='UTF-8') as f: + guan.write_two_dimensional_data_without_xy_array_and_without_opening_file(matrix, f) + guan.statistics_of_guan_package() + +# 在文件中写入二维数据（不包括x和y）（需要输入文件） +def write_two_dimensional_data_without_xy_array_and_without_opening_file(matrix, f): + for row in matrix: + for element in row: + f.write(str(element)+' ') + f.write('\n') + import guan + guan.statistics_of_guan_package() + +# 以显示编号的样式，打印数组 +def print_array_with_index(array, show_index=1, index_type=0): + if show_index==0: + for i0 in array: + print(i0) + else: + if index_type==0: + index = 0 + for i0 in array: + print(index, i0) + index += 1 + else: + index = 0 + for i0 in array: + index += 1 + print(index, i0) + import guan + guan.statistics_of_guan_package() diff --git a/PyPI/src/guan/topological_invariant.py b/PyPI/src/guan/topological_invariant.py new file mode 100644 index 0000000..325c5c6 --- /dev/null +++ b/PyPI/src/guan/topological_invariant.py @@ -0,0 +1,552 @@ +# Module: topological_invariant + +# 通过高效法计算方格子的陈数 +def calculate_chern_number_for_square_lattice_with_efficient_method(hamiltonian_function, precision=100, print_show=0): + import numpy as np + import math + import cmath + import guan + if np.array(hamiltonian_function(0, 0)).shape==(): + dim = 1 + else: + dim = np.array(hamiltonian_function(0, 0)).shape[0] + delta = 2*math.pi/precision + chern_number = np.zeros(dim, dtype=complex) + for kx in np.arange(-math.pi, math.pi, delta): + if print_show == 1: + print(kx) + for ky in np.arange(-math.pi, math.pi, delta): + H = hamiltonian_function(kx, ky) + vector = guan.calculate_eigenvector(H) + H_delta_kx = hamiltonian_function(kx+delta, ky) + vector_delta_kx = guan.calculate_eigenvector(H_delta_kx) + H_delta_ky = hamiltonian_function(kx, ky+delta) + vector_delta_ky = guan.calculate_eigenvector(H_delta_ky) + H_delta_kx_ky = hamiltonian_function(kx+delta, ky+delta) + vector_delta_kx_ky = guan.calculate_eigenvector(H_delta_kx_ky) + for i in range(dim): + vector_i = vector[:, i] + vector_delta_kx_i = vector_delta_kx[:, i] + vector_delta_ky_i = vector_delta_ky[:, i] + vector_delta_kx_ky_i = vector_delta_kx_ky[:, i] + Ux = np.dot(np.conj(vector_i), vector_delta_kx_i)/abs(np.dot(np.conj(vector_i), vector_delta_kx_i)) + Uy = np.dot(np.conj(vector_i), vector_delta_ky_i)/abs(np.dot(np.conj(vector_i), vector_delta_ky_i)) + Ux_y = np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i)) + Uy_x = np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i)) + F = cmath.log(Ux*Uy_x*(1/Ux_y)*(1/Uy)) + chern_number[i] = chern_number[i] + F + chern_number = chern_number/(2*math.pi*1j) + guan.statistics_of_guan_package() + return chern_number + +# 通过高效法计算方格子的陈数（可计算简并的情况） +def calculate_chern_number_for_square_lattice_with_efficient_method_for_degenerate_case(hamiltonian_function, index_of_bands=[0, 1], precision=100, print_show=0): + import numpy as np + import math + import cmath + delta = 2*math.pi/precision + chern_number = 0 + for kx in np.arange(-math.pi, math.pi, delta): + if print_show == 1: + print(kx) + for ky in np.arange(-math.pi, math.pi, delta): + H = hamiltonian_function(kx, ky) + eigenvalue, vector = np.linalg.eigh(H) + H_delta_kx = hamiltonian_function(kx+delta, ky) + eigenvalue, vector_delta_kx = np.linalg.eigh(H_delta_kx) + H_delta_ky = hamiltonian_function(kx, ky+delta) + eigenvalue, vector_delta_ky = np.linalg.eigh(H_delta_ky) + H_delta_kx_ky = hamiltonian_function(kx+delta, ky+delta) + eigenvalue, vector_delta_kx_ky = np.linalg.eigh(H_delta_kx_ky) + dim = len(index_of_bands) + det_value = 1 + # first dot product + dot_matrix = np.zeros((dim , dim), dtype=complex) + i0 = 0 + for dim1 in index_of_bands: + j0 = 0 + for dim2 in index_of_bands: + dot_matrix[i0, j0] = np.dot(np.conj(vector[:, dim1]), vector_delta_kx[:, dim2]) + j0 += 1 + i0 += 1 + dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) + det_value = det_value*dot_matrix + # second dot product + dot_matrix = np.zeros((dim , dim), dtype=complex) + i0 = 0 + for dim1 in index_of_bands: + j0 = 0 + for dim2 in index_of_bands: + dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx[:, dim1]), vector_delta_kx_ky[:, dim2]) + j0 += 1 + i0 += 1 + dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) + det_value = det_value*dot_matrix + # third dot product + dot_matrix = np.zeros((dim , dim), dtype=complex) + i0 = 0 + for dim1 in index_of_bands: + j0 = 0 + for dim2 in index_of_bands: + dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx_ky[:, dim1]), vector_delta_ky[:, dim2]) + j0 += 1 + i0 += 1 + dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) + det_value = det_value*dot_matrix + # four dot product + dot_matrix = np.zeros((dim , dim), dtype=complex) + i0 = 0 + for dim1 in index_of_bands: + j0 = 0 + for dim2 in index_of_bands: + dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_ky[:, dim1]), vector[:, dim2]) + j0 += 1 + i0 += 1 + dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) + det_value= det_value*dot_matrix + chern_number += cmath.log(det_value) + chern_number = chern_number/(2*math.pi*1j) + import guan + guan.statistics_of_guan_package() + return chern_number + +# 通过Wilson loop方法计算方格子的陈数 +def calculate_chern_number_for_square_lattice_with_wilson_loop(hamiltonian_function, precision_of_plaquettes=20, precision_of_wilson_loop=5, print_show=0): + import numpy as np + import math + delta = 2*math.pi/precision_of_plaquettes + chern_number = 0 + for kx in np.arange(-math.pi, math.pi, delta): + if print_show == 1: + print(kx) + for ky in np.arange(-math.pi, math.pi, delta): + vector_array = [] + # line_1 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx+delta/precision_of_wilson_loop*i0, ky) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + # line_2 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx+delta, ky+delta/precision_of_wilson_loop*i0) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + # line_3 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx+delta-delta/precision_of_wilson_loop*i0, ky+delta) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + # line_4 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx, ky+delta-delta/precision_of_wilson_loop*i0) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + wilson_loop = 1 + for i0 in range(len(vector_array)-1): + wilson_loop = wilson_loop*np.dot(vector_array[i0].transpose().conj(), vector_array[i0+1]) + wilson_loop = wilson_loop*np.dot(vector_array[len(vector_array)-1].transpose().conj(), vector_array[0]) + arg = np.log(np.diagonal(wilson_loop))/1j + chern_number = chern_number + arg + chern_number = chern_number/(2*math.pi) + import guan + guan.statistics_of_guan_package() + return chern_number + +# 通过Wilson loop方法计算方格子的陈数（可计算简并的情况） +def calculate_chern_number_for_square_lattice_with_wilson_loop_for_degenerate_case(hamiltonian_function, index_of_bands=[0, 1], precision_of_plaquettes=20, precision_of_wilson_loop=5, print_show=0): + import numpy as np + import math + delta = 2*math.pi/precision_of_plaquettes + chern_number = 0 + for kx in np.arange(-math.pi, math.pi, delta): + if print_show == 1: + print(kx) + for ky in np.arange(-math.pi, math.pi, delta): + vector_array = [] + # line_1 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx+delta/precision_of_wilson_loop*i0, ky) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + # line_2 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx+delta, ky+delta/precision_of_wilson_loop*i0) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + # line_3 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx+delta-delta/precision_of_wilson_loop*i0, ky+delta) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + # line_4 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx, ky+delta-delta/precision_of_wilson_loop*i0) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + wilson_loop = 1 + dim = len(index_of_bands) + for i0 in range(len(vector_array)-1): + dot_matrix = np.zeros((dim , dim), dtype=complex) + i01 = 0 + for dim1 in index_of_bands: + i02 = 0 + for dim2 in index_of_bands: + dot_matrix[i01, i02] = np.dot(vector_array[i0][:, dim1].transpose().conj(), vector_array[i0+1][:, dim2]) + i02 += 1 + i01 += 1 + det_value = np.linalg.det(dot_matrix) + wilson_loop = wilson_loop*det_value + dot_matrix_plus = np.zeros((dim , dim), dtype=complex) + i01 = 0 + for dim1 in index_of_bands: + i02 = 0 + for dim2 in index_of_bands: + dot_matrix_plus[i01, i02] = np.dot(vector_array[len(vector_array)-1][:, dim1].transpose().conj(), vector_array[0][:, dim2]) + i02 += 1 + i01 += 1 + det_value = np.linalg.det(dot_matrix_plus) + wilson_loop = wilson_loop*det_value + arg = np.log(wilson_loop)/1j + chern_number = chern_number + arg + chern_number = chern_number/(2*math.pi) + import guan + guan.statistics_of_guan_package() + return chern_number + +# 通过高效法计算贝利曲率 +def calculate_berry_curvature_with_efficient_method(hamiltonian_function, k_min='default', k_max='default', precision=100, print_show=0): + import numpy as np + import cmath + import guan + import math + if k_min == 'default': + k_min = -math.pi + if k_max == 'default': + k_max=math.pi + if np.array(hamiltonian_function(0, 0)).shape==(): + dim = 1 + else: + dim = np.array(hamiltonian_function(0, 0)).shape[0] + delta = (k_max-k_min)/precision + k_array = np.arange(k_min, k_max, delta) + berry_curvature_array = np.zeros((k_array.shape[0], k_array.shape[0], dim), dtype=complex) + i0 = 0 + for kx in k_array: + if print_show == 1: + print(kx) + j0 = 0 + for ky in k_array: + H = hamiltonian_function(kx, ky) + vector = guan.calculate_eigenvector(H) + H_delta_kx = hamiltonian_function(kx+delta, ky) + vector_delta_kx = guan.calculate_eigenvector(H_delta_kx) + H_delta_ky = hamiltonian_function(kx, ky+delta) + vector_delta_ky = guan.calculate_eigenvector(H_delta_ky) + H_delta_kx_ky = hamiltonian_function(kx+delta, ky+delta) + vector_delta_kx_ky = guan.calculate_eigenvector(H_delta_kx_ky) + for i in range(dim): + vector_i = vector[:, i] + vector_delta_kx_i = vector_delta_kx[:, i] + vector_delta_ky_i = vector_delta_ky[:, i] + vector_delta_kx_ky_i = vector_delta_kx_ky[:, i] + Ux = np.dot(np.conj(vector_i), vector_delta_kx_i)/abs(np.dot(np.conj(vector_i), vector_delta_kx_i)) + Uy = np.dot(np.conj(vector_i), vector_delta_ky_i)/abs(np.dot(np.conj(vector_i), vector_delta_ky_i)) + Ux_y = np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i)) + Uy_x = np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i)) + berry_curvature = cmath.log(Ux*Uy_x*(1/Ux_y)*(1/Uy))/delta/delta*1j + berry_curvature_array[j0, i0, i] = berry_curvature + j0 += 1 + i0 += 1 + guan.statistics_of_guan_package() + return k_array, berry_curvature_array + +# 通过高效法计算贝利曲率（可计算简并的情况） +def calculate_berry_curvature_with_efficient_method_for_degenerate_case(hamiltonian_function, index_of_bands=[0, 1], k_min='default', k_max='default', precision=100, print_show=0): + import numpy as np + import cmath + import math + if k_min == 'default': + k_min = -math.pi + if k_max == 'default': + k_max=math.pi + delta = (k_max-k_min)/precision + k_array = np.arange(k_min, k_max, delta) + berry_curvature_array = np.zeros((k_array.shape[0], k_array.shape[0]), dtype=complex) + i00 = 0 + for kx in np.arange(k_min, k_max, delta): + if print_show == 1: + print(kx) + j00 = 0 + for ky in np.arange(k_min, k_max, delta): + H = hamiltonian_function(kx, ky) + eigenvalue, vector = np.linalg.eigh(H) + H_delta_kx = hamiltonian_function(kx+delta, ky) + eigenvalue, vector_delta_kx = np.linalg.eigh(H_delta_kx) + H_delta_ky = hamiltonian_function(kx, ky+delta) + eigenvalue, vector_delta_ky = np.linalg.eigh(H_delta_ky) + H_delta_kx_ky = hamiltonian_function(kx+delta, ky+delta) + eigenvalue, vector_delta_kx_ky = np.linalg.eigh(H_delta_kx_ky) + dim = len(index_of_bands) + det_value = 1 + # first dot product + dot_matrix = np.zeros((dim , dim), dtype=complex) + i0 = 0 + for dim1 in index_of_bands: + j0 = 0 + for dim2 in index_of_bands: + dot_matrix[i0, j0] = np.dot(np.conj(vector[:, dim1]), vector_delta_kx[:, dim2]) + j0 += 1 + i0 += 1 + dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) + det_value = det_value*dot_matrix + # second dot product + dot_matrix = np.zeros((dim , dim), dtype=complex) + i0 = 0 + for dim1 in index_of_bands: + j0 = 0 + for dim2 in index_of_bands: + dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx[:, dim1]), vector_delta_kx_ky[:, dim2]) + j0 += 1 + i0 += 1 + dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) + det_value = det_value*dot_matrix + # third dot product + dot_matrix = np.zeros((dim , dim), dtype=complex) + i0 = 0 + for dim1 in index_of_bands: + j0 = 0 + for dim2 in index_of_bands: + dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx_ky[:, dim1]), vector_delta_ky[:, dim2]) + j0 += 1 + i0 += 1 + dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) + det_value = det_value*dot_matrix + # four dot product + dot_matrix = np.zeros((dim , dim), dtype=complex) + i0 = 0 + for dim1 in index_of_bands: + j0 = 0 + for dim2 in index_of_bands: + dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_ky[:, dim1]), vector[:, dim2]) + j0 += 1 + i0 += 1 + dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix)) + det_value= det_value*dot_matrix + berry_curvature = cmath.log(det_value)/delta/delta*1j + berry_curvature_array[j00, i00] = berry_curvature + j00 += 1 + i00 += 1 + import guan + guan.statistics_of_guan_package() + return k_array, berry_curvature_array + +# 通过Wilson loop方法计算贝里曲率 +def calculate_berry_curvature_with_wilson_loop(hamiltonian_function, k_min='default', k_max='default', precision_of_plaquettes=20, precision_of_wilson_loop=5, print_show=0): + import numpy as np + import math + if k_min == 'default': + k_min = -math.pi + if k_max == 'default': + k_max=math.pi + if np.array(hamiltonian_function(0, 0)).shape==(): + dim = 1 + else: + dim = np.array(hamiltonian_function(0, 0)).shape[0] + delta = (k_max-k_min)/precision_of_plaquettes + k_array = np.arange(k_min, k_max, delta) + berry_curvature_array = np.zeros((k_array.shape[0], k_array.shape[0], dim), dtype=complex) + i00 = 0 + for kx in k_array: + if print_show == 1: + print(kx) + j00 = 0 + for ky in k_array: + vector_array = [] + # line_1 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx+delta/precision_of_wilson_loop*i0, ky) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + # line_2 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx+delta, ky+delta/precision_of_wilson_loop*i0) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + # line_3 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx+delta-delta/precision_of_wilson_loop*i0, ky+delta) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + # line_4 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx, ky+delta-delta/precision_of_wilson_loop*i0) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + wilson_loop = 1 + for i0 in range(len(vector_array)-1): + wilson_loop = wilson_loop*np.dot(vector_array[i0].transpose().conj(), vector_array[i0+1]) + wilson_loop = wilson_loop*np.dot(vector_array[len(vector_array)-1].transpose().conj(), vector_array[0]) + berry_curvature = np.log(np.diagonal(wilson_loop))/delta/delta*1j + berry_curvature_array[j00, i00, :]=berry_curvature + j00 += 1 + i00 += 1 + import guan + guan.statistics_of_guan_package() + return k_array, berry_curvature_array + +# 通过Wilson loop方法计算贝里曲率（可计算简并的情况） +def calculate_berry_curvature_with_wilson_loop_for_degenerate_case(hamiltonian_function, index_of_bands=[0, 1], k_min='default', k_max='default', precision_of_plaquettes=20, precision_of_wilson_loop=5, print_show=0): + import numpy as np + import math + if k_min == 'default': + k_min = -math.pi + if k_max == 'default': + k_max=math.pi + delta = (k_max-k_min)/precision_of_plaquettes + k_array = np.arange(k_min, k_max, delta) + berry_curvature_array = np.zeros((k_array.shape[0], k_array.shape[0]), dtype=complex) + i000 = 0 + for kx in k_array: + if print_show == 1: + print(kx) + j000 = 0 + for ky in k_array: + vector_array = [] + # line_1 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx+delta/precision_of_wilson_loop*i0, ky) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + # line_2 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx+delta, ky+delta/precision_of_wilson_loop*i0) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + # line_3 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx+delta-delta/precision_of_wilson_loop*i0, ky+delta) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + # line_4 + for i0 in range(precision_of_wilson_loop): + H_delta = hamiltonian_function(kx, ky+delta-delta/precision_of_wilson_loop*i0) + eigenvalue, eigenvector = np.linalg.eig(H_delta) + vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))] + vector_array.append(vector_delta) + wilson_loop = 1 + dim = len(index_of_bands) + for i0 in range(len(vector_array)-1): + dot_matrix = np.zeros((dim , dim), dtype=complex) + i01 = 0 + for dim1 in index_of_bands: + i02 = 0 + for dim2 in index_of_bands: + dot_matrix[i01, i02] = np.dot(vector_array[i0][:, dim1].transpose().conj(), vector_array[i0+1][:, dim2]) + i02 += 1 + i01 += 1 + det_value = np.linalg.det(dot_matrix) + wilson_loop = wilson_loop*det_value + dot_matrix_plus = np.zeros((dim , dim), dtype=complex) + i01 = 0 + for dim1 in index_of_bands: + i02 = 0 + for dim2 in index_of_bands: + dot_matrix_plus[i01, i02] = np.dot(vector_array[len(vector_array)-1][:, dim1].transpose().conj(), vector_array[0][:, dim2]) + i02 += 1 + i01 += 1 + det_value = np.linalg.det(dot_matrix_plus) + wilson_loop = wilson_loop*det_value + berry_curvature = np.log(wilson_loop)/delta/delta*1j + berry_curvature_array[j000, i000]=berry_curvature + j000 += 1 + i000 += 1 + import guan + guan.statistics_of_guan_package() + return k_array, berry_curvature_array + +# 计算蜂窝格子的陈数（高效法） +def calculate_chern_number_for_honeycomb_lattice(hamiltonian_function, a=1, precision=300, print_show=0): + import numpy as np + import math + import cmath + import guan + if np.array(hamiltonian_function(0, 0)).shape==(): + dim = 1 + else: + dim = np.array(hamiltonian_function(0, 0)).shape[0] + chern_number = np.zeros(dim, dtype=complex) + L1 = 4*math.sqrt(3)*math.pi/9/a + L2 = 2*math.sqrt(3)*math.pi/9/a + L3 = 2*math.pi/3/a + delta1 = 2*L1/precision + delta3 = 2*L3/precision + for kx in np.arange(-L1, L1, delta1): + if print_show == 1: + print(kx) + for ky in np.arange(-L3, L3, delta3): + if (-L2<=kx<=L2) or (kx>L2 and -(L1-kx)*math.tan(math.pi/3)<=ky<=(L1-kx)*math.tan(math.pi/3)) or (kx<-L2 and -(kx-(-L1))*math.tan(math.pi/3)<=ky<=(kx-(-L1))*math.tan(math.pi/3)): + H = hamiltonian_function(kx, ky) + vector = guan.calculate_eigenvector(H) + H_delta_kx = hamiltonian_function(kx+delta1, ky) + vector_delta_kx = guan.calculate_eigenvector(H_delta_kx) + H_delta_ky = hamiltonian_function(kx, ky+delta3) + vector_delta_ky = guan.calculate_eigenvector(H_delta_ky) + H_delta_kx_ky = hamiltonian_function(kx+delta1, ky+delta3) + vector_delta_kx_ky = guan.calculate_eigenvector(H_delta_kx_ky) + for i in range(dim): + vector_i = vector[:, i] + vector_delta_kx_i = vector_delta_kx[:, i] + vector_delta_ky_i = vector_delta_ky[:, i] + vector_delta_kx_ky_i = vector_delta_kx_ky[:, i] + Ux = np.dot(np.conj(vector_i), vector_delta_kx_i)/abs(np.dot(np.conj(vector_i), vector_delta_kx_i)) + Uy = np.dot(np.conj(vector_i), vector_delta_ky_i)/abs(np.dot(np.conj(vector_i), vector_delta_ky_i)) + Ux_y = np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i)) + Uy_x = np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i)) + F = cmath.log(Ux*Uy_x*(1/Ux_y)*(1/Uy)) + chern_number[i] = chern_number[i] + F + chern_number = chern_number/(2*math.pi*1j) + guan.statistics_of_guan_package() + return chern_number + +# 计算Wilson loop +def calculate_wilson_loop(hamiltonian_function, k_min='default', k_max='default', precision=100, print_show=0): + import numpy as np + import guan + import math + if k_min == 'default': + k_min = -math.pi + if k_max == 'default': + k_max=math.pi + k_array = np.linspace(k_min, k_max, precision) + dim = np.array(hamiltonian_function(0)).shape[0] + wilson_loop_array = np.ones(dim, dtype=complex) + for i in range(dim): + if print_show == 1: + print(i) + eigenvector_array = [] + for k in k_array: + eigenvector = guan.calculate_eigenvector(hamiltonian_function(k)) + if k != k_max: + eigenvector_array.append(eigenvector[:, i]) + else: + eigenvector_array.append(eigenvector_array[0]) + for i0 in range(precision-1): + F = np.dot(eigenvector_array[i0+1].transpose().conj(), eigenvector_array[i0]) + wilson_loop_array[i] = np.dot(F, wilson_loop_array[i]) + guan.statistics_of_guan_package() + return wilson_loop_array \ No newline at end of file diff --git a/README.md b/README.md index 7306f29..9c0ce0a 100755 --- a/README.md +++ b/README.md @@ -17,7 +17,7 @@ import guan + basic functions + Fourier transform + Hamiltonian of finite size systems -+ Hamiltonian of models in the reciprocal space ++ Hamiltonian of models in reciprocal space + band structures and wave functions + Green functions + density of states