diff --git a/API_Reference/API_Reference.py b/API_Reference/API_Reference.py index df64c89..f03f58a 100644 --- a/API_Reference/API_Reference.py +++ b/API_Reference/API_Reference.py @@ -27,6 +27,9 @@ sigma_zz = guan.sigma_zz() hamiltonian = guan.one_dimensional_fourier_transform(k, unit_cell, hopping) hamiltonian = guan.two_dimensional_fourier_transform_for_square_lattice(k1, k2, unit_cell, hopping_1, hopping_2) hamiltonian = guan.three_dimensional_fourier_transform_for_cubic_lattice(k1, k2, k3, unit_cell, hopping_1, hopping_2, hopping_3) +hamiltonian_function = guan.one_dimensional_fourier_transform_with_k(unit_cell, hopping) +hamiltonian_function = guan.two_dimensional_fourier_transform_for_square_lattice_with_k1_k2(unit_cell, hopping_1, hopping_2) +hamiltonian_function = guan.three_dimensional_fourier_transform_for_cubic_lattice_with_k1_k2_k3(unit_cell, hopping_1, hopping_2, hopping_3) b1 = guan.calculate_one_dimensional_reciprocal_lattice_vector(a1) b1, b2 = guan.calculate_two_dimensional_reciprocal_lattice_vectors(a1, a2) b1, b2, b3 = guan.calculate_three_dimensional_reciprocal_lattice_vectors(a1, a2, a3) @@ -35,11 +38,11 @@ b1, b2 = guan.calculate_two_dimensional_reciprocal_lattice_vectors_with_sympy(a1 b1, b2, b3 = guan.calculate_three_dimensional_reciprocal_lattice_vectors_with_sympy(a1, a2, a3) # Hamiltonian of finite size systems -hamiltonian = guan.finite_size_along_one_direction(N, on_site=0, hopping=1, period=0) -hamiltonian = guan.finite_size_along_two_directions_for_square_lattice(N1, N2, on_site=0, hopping_1=1, hopping_2=1, period_1=0, period_2=0) -hamiltonian = guan.finite_size_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) -hopping = guan.hopping_along_zigzag_direction_for_graphene(N) -hamiltonian = guan.finite_size_along_two_directions_for_graphene(N1, N2, period_1=0, period_2=0) +hamiltonian = guan.hamiltonian_of_finite_size_system_along_one_direction(N, on_site=0, hopping=1, period=0) +hamiltonian = guan.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) +hamiltonian = guan.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) +hopping = guan.hopping_matrix_along_zigzag_direction_for_graphene_ribbon(N) +hamiltonian = guan.hamiltonian_of_finite_size_system_along_two_directions_for_graphene(N1, N2, period_1=0, period_2=0) # Hamiltonian of models in the reciprocal space hamiltonian = guan.hamiltonian_of_simple_chain(k) @@ -51,6 +54,7 @@ hamiltonian = guan.hamiltonian_of_graphene(k1, k2, M=0, t=1, a=1/sqrt(3)) hamiltonian = guan.hamiltonian_of_graphene_with_zigzag_in_quasi_one_dimension(k, N=10, M=0, t=1) hamiltonian = guan.hamiltonian_of_haldane_model(k1, k2, M=2/3, t1=1, t2=1/3, phi=pi/4, a=1/sqrt(3)) hamiltonian = guan.hamiltonian_of_haldane_model_in_quasi_one_dimension(k, N=10, M=2/3, t1=1, t2=1/3, phi=pi/4) +hamiltonian = guan.hamiltonian_of_one_QAH_model(k1, k2, t1=1, t2=1, t3=0.5, m=-1) # band structures and wave functions eigenvalue = guan.calculate_eigenvalue(hamiltonian) diff --git a/PyPI/setup.cfg b/PyPI/setup.cfg index f5c7aa7..63f44bc 100644 --- a/PyPI/setup.cfg +++ b/PyPI/setup.cfg @@ -1,7 +1,7 @@ [metadata] # replace with your username: name = guan -version = 0.0.49 +version = 0.0.51 author = guanjihuan author_email = guanjihuan@163.com description = An open source python package diff --git a/PyPI/src/guan/Fourier_transform.py b/PyPI/src/guan/Fourier_transform.py deleted file mode 100644 index cd52519..0000000 --- a/PyPI/src/guan/Fourier_transform.py +++ /dev/null @@ -1,85 +0,0 @@ -# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about. The primary location of this package is on website https://py.guanjihuan.com. - -# Fourier_transform - -import numpy as np -import cmath - -# Fourier_transform for discrete lattices - -def one_dimensional_fourier_transform(k, unit_cell, hopping): - 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) - return hamiltonian - -def two_dimensional_fourier_transform_for_square_lattice(k1, k2, unit_cell, hopping_1, hopping_2): - 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) - return hamiltonian - -def three_dimensional_fourier_transform_for_cubic_lattice(k1, k2, k3, unit_cell, hopping_1, hopping_2, hopping_3): - 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) - return hamiltonian - - -## calculate reciprocal lattice vectors - -def calculate_one_dimensional_reciprocal_lattice_vector(a1): - b1 = 2*np.pi/a1 - return b1 - -def calculate_two_dimensional_reciprocal_lattice_vectors(a1, a2): - 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) - return b1, b2 - -def calculate_three_dimensional_reciprocal_lattice_vectors(a1, a2, a3): - 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)) - return b1, b2, b3 - -def calculate_one_dimensional_reciprocal_lattice_vector_with_sympy(a1): - import sympy - b1 = 2*sympy.pi/a1 - 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]]) - 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) - return b1, b2, b3 \ No newline at end of file diff --git a/PyPI/src/guan/Green_functions.py b/PyPI/src/guan/Green_functions.py deleted file mode 100644 index 17dbbf8..0000000 --- a/PyPI/src/guan/Green_functions.py +++ /dev/null @@ -1,110 +0,0 @@ -# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about. The primary location of this package is on website https://py.guanjihuan.com. - -# Green functions - -import numpy as np -import guan - -def green_function(fermi_energy, hamiltonian, broadening, self_energy=0): - 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) - return green - -def green_function_nn_n(fermi_energy, h00, h01, green_nn_n_minus, broadening, self_energy=0): - 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) - return green_nn_n - -def green_function_in_n(green_in_n_minus, h01, green_nn_n): - green_in_n = np.dot(np.dot(green_in_n_minus, h01), green_nn_n) - return green_in_n - -def green_function_ni_n(green_nn_n, h01, green_ni_n_minus): - h01 = np.array(h01) - green_ni_n = np.dot(np.dot(green_nn_n, h01.transpose().conj()), green_ni_n_minus) - return green_ni_n - -def green_function_ii_n(green_ii_n_minus, green_in_n_minus, h01, green_nn_n, green_ni_n_minus): - 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) - return green_ii_n - -def transfer_matrix(fermi_energy, h00, h01): - 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 - return transfer - -def surface_green_function_of_lead(fermi_energy, h00, h01): - 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))) - return right_lead_surface, left_lead_surface - -def self_energy_of_lead(fermi_energy, h00, h01): - 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 - 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): - 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 - 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): - 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 - return self_energy, gamma - -def green_function_with_leads(fermi_energy, h00, h01, h_LC, h_CR, center_hamiltonian): - 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) - return green, gamma_right, gamma_left - -def electron_correlation_function_green_n_for_local_current(fermi_energy, h00, h01, h_LC, h_CR, center_hamiltonian): - 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())) - 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 deleted file mode 100644 index 17e820c..0000000 --- a/PyPI/src/guan/Hamiltonian_of_finite_size_systems.py +++ /dev/null @@ -1,115 +0,0 @@ -# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about. The primary location of this package is on website https://py.guanjihuan.com. - -# Hamiltonian of finite size systems - -import numpy as np -import guan - -def finite_size_along_one_direction(N, on_site=0, hopping=1, period=0): - 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() - return hamiltonian - -def finite_size_along_two_directions_for_square_lattice(N1, N2, on_site=0, hopping_1=1, hopping_2=1, period_1=0, period_2=0): - 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() - return hamiltonian - -def finite_size_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): - 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() - return hamiltonian - -def hopping_along_zigzag_direction_for_graphene(N): - hopping = np.zeros((4*N, 4*N), dtype=complex) - for i0 in range(N): - hopping[4*i0+1, 4*i0+0] = 1 - hopping[4*i0+2, 4*i0+3] = 1 - return hopping - -def finite_size_along_two_directions_for_graphene(N1, N2, period_1=0, period_2=0): - on_site = guan.finite_size_along_one_direction(4) - hopping_1 = guan.hopping_along_zigzag_direction_for_graphene(1) - hopping_2 = np.zeros((4, 4), dtype=complex) - hopping_2[3, 0] = 1 - hamiltonian = guan.finite_size_along_two_directions_for_square_lattice(N1, N2, on_site, hopping_1, hopping_2, period_1, period_2) - return hamiltonian \ No newline at end of file diff --git a/PyPI/src/guan/Hamiltonian_of_models_in_the_reciprocal_space.py b/PyPI/src/guan/Hamiltonian_of_models_in_the_reciprocal_space.py deleted file mode 100644 index ff2ce0d..0000000 --- a/PyPI/src/guan/Hamiltonian_of_models_in_the_reciprocal_space.py +++ /dev/null @@ -1,124 +0,0 @@ -# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about. The primary location of this package is on website https://py.guanjihuan.com. - -# Hamiltonian of models in the reciprocal space - -import numpy as np -import cmath -from math import * -import guan - -def hamiltonian_of_simple_chain(k): - hamiltonian = guan.one_dimensional_fourier_transform(k, unit_cell=0, hopping=1) - return hamiltonian - -def hamiltonian_of_square_lattice(k1, k2): - hamiltonian = guan.two_dimensional_fourier_transform_for_square_lattice(k1, k2, unit_cell=0, hopping_1=1, hopping_2=1) - return hamiltonian - -def hamiltonian_of_square_lattice_in_quasi_one_dimension(k, N=10): - 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 - for i in range(N): - h01[i, i] = 1 - hamiltonian = guan.one_dimensional_fourier_transform(k, unit_cell=h00, hopping=h01) - return hamiltonian - -def hamiltonian_of_cubic_lattice(k1, k2, k3): - hamiltonian = guan.three_dimensional_fourier_transform_for_cubic_lattice(k1, k2, k3, unit_cell=0, hopping_1=1, hopping_2=1, hopping_3=1) - return hamiltonian - -def hamiltonian_of_ssh_model(k, v=0.6, w=1): - 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) - return hamiltonian - -def hamiltonian_of_graphene(k1, k2, M=0, t=1, a=1/sqrt(3)): - h0 = np.zeros((2, 2), dtype=complex) # mass term - h1 = np.zeros((2, 2), dtype=complex) # nearest hopping - h0[0, 0] = M - h0[1, 1] = -M - h1[1, 0] = t*(cmath.exp(1j*k2*a)+cmath.exp(1j*sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j/2*k2*a)) - h1[0, 1] = h1[1, 0].conj() - hamiltonian = h0 + h1 - return hamiltonian - -def hamiltonian_of_graphene_with_zigzag_in_quasi_one_dimension(k, N=10, M=0, t=1): - 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 - 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) - return hamiltonian - -def hamiltonian_of_haldane_model(k1, k2, M=2/3, t1=1, t2=1/3, phi=pi/4, a=1/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*sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*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*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a)) - h2[1, 1] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a)) - hamiltonian = h0 + h1 + h2 + h2.transpose().conj() - return hamiltonian - -def hamiltonian_of_haldane_model_in_quasi_one_dimension(k, N=10, M=2/3, t1=1, t2=1/3, phi=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() - 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) - return hamiltonian \ No newline at end of file diff --git a/PyPI/src/guan/__init__.py b/PyPI/src/guan/__init__.py index 0321edd..ea43d13 100644 --- a/PyPI/src/guan/__init__.py +++ b/PyPI/src/guan/__init__.py @@ -1,16 +1,1568 @@ # Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about. The primary location of this package is on website https://py.guanjihuan.com. -# import modules +# Modules -from .basic_functions import * -from .Fourier_transform import * -from .Hamiltonian_of_finite_size_systems import * -from .Hamiltonian_of_models_in_the_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 .read_and_write import * -from .plot_figures import * -from .others import * \ No newline at end of file +# # 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: read_and_write +# # Module 11: plot_figures +# # Module 12: others + +# import packages + +import numpy as np +from math import * +import cmath +import functools +import copy +import guan + +# Module 1: basic functions + +## test + +def test(): + print('\nSuccess in the installation of Guan package!\n') + +## Pauli matrices + +def sigma_0(): + return np.eye(2) + +def sigma_x(): + return np.array([[0, 1],[1, 0]]) + +def sigma_y(): + return np.array([[0, -1j],[1j, 0]]) + +def sigma_z(): + return np.array([[1, 0],[0, -1]]) + +## Kronecker product of Pauli matrices + +def sigma_00(): + return np.kron(sigma_0(), sigma_0()) + +def sigma_0x(): + return np.kron(sigma_0(), sigma_x()) + +def sigma_0y(): + return np.kron(sigma_0(), sigma_y()) + +def sigma_0z(): + return np.kron(sigma_0(), sigma_z()) + +def sigma_x0(): + return np.kron(sigma_x(), sigma_0()) + +def sigma_xx(): + return np.kron(sigma_x(), sigma_x()) + +def sigma_xy(): + return np.kron(sigma_x(), sigma_y()) + +def sigma_xz(): + return np.kron(sigma_x(), sigma_z()) + +def sigma_y0(): + return np.kron(sigma_y(), sigma_0()) + +def sigma_yx(): + return np.kron(sigma_y(), sigma_x()) + +def sigma_yy(): + return np.kron(sigma_y(), sigma_y()) + +def sigma_yz(): + return np.kron(sigma_y(), sigma_z()) + +def sigma_z0(): + return np.kron(sigma_z(), sigma_0()) + +def sigma_zx(): + return np.kron(sigma_z(), sigma_x()) + +def sigma_zy(): + return np.kron(sigma_z(), sigma_y()) + +def sigma_zz(): + return np.kron(sigma_z(), sigma_z()) + + + + + + + + + + +# Module 2: Fourier_transform + +# Fourier_transform for discrete lattices + +def one_dimensional_fourier_transform(k, unit_cell, hopping): + 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) + return hamiltonian + +def two_dimensional_fourier_transform_for_square_lattice(k1, k2, unit_cell, hopping_1, hopping_2): + 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) + return hamiltonian + +def three_dimensional_fourier_transform_for_cubic_lattice(k1, k2, k3, unit_cell, hopping_1, hopping_2, hopping_3): + 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) + return hamiltonian + +def one_dimensional_fourier_transform_with_k(unit_cell, hopping): + hamiltonian_function = functools.partial(guan.one_dimensional_fourier_transform, unit_cell=unit_cell, hopping=hopping) + return hamiltonian_function + +def two_dimensional_fourier_transform_for_square_lattice_with_k1_k2(unit_cell, hopping_1, hopping_2): + hamiltonian_function = functools.partial(guan.two_dimensional_fourier_transform_for_square_lattice, unit_cell=unit_cell, hopping_1=hopping_1, hopping_2=hopping_2) + return hamiltonian_function + +def three_dimensional_fourier_transform_for_cubic_lattice_with_k1_k2_k3(unit_cell, hopping_1, hopping_2, hopping_3): + 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) + return hamiltonian_function + +## calculate reciprocal lattice vectors + +def calculate_one_dimensional_reciprocal_lattice_vector(a1): + b1 = 2*np.pi/a1 + return b1 + +def calculate_two_dimensional_reciprocal_lattice_vectors(a1, a2): + 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) + return b1, b2 + +def calculate_three_dimensional_reciprocal_lattice_vectors(a1, a2, a3): + 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)) + return b1, b2, b3 + +def calculate_one_dimensional_reciprocal_lattice_vector_with_sympy(a1): + import sympy + b1 = 2*sympy.pi/a1 + 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]]) + 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) + 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): + 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() + 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): + 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() + 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): + 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() + return hamiltonian + +def hopping_matrix_along_zigzag_direction_for_graphene_ribbon(N): + hopping = np.zeros((4*N, 4*N), dtype=complex) + for i0 in range(N): + hopping[4*i0+1, 4*i0+0] = 1 + hopping[4*i0+2, 4*i0+3] = 1 + return hopping + +def hamiltonian_of_finite_size_system_along_two_directions_for_graphene(N1, N2, period_1=0, period_2=0): + on_site = guan.hamiltonian_of_finite_size_system_along_one_direction(4) + hopping_1 = guan.hopping_matrix_along_zigzag_direction_for_graphene_ribbon(1) + hopping_2 = np.zeros((4, 4), dtype=complex) + hopping_2[3, 0] = 1 + hamiltonian = guan.finite_size_along_two_directions_for_square_lattice(N1, N2, on_site, hopping_1, hopping_2, period_1, period_2) + return hamiltonian + + + + + + + + + + +# Module 4: Hamiltonian of models in the reciprocal space + +def hamiltonian_of_simple_chain(k): + hamiltonian = guan.one_dimensional_fourier_transform(k, unit_cell=0, hopping=1) + return hamiltonian + +def hamiltonian_of_square_lattice(k1, k2): + hamiltonian = guan.two_dimensional_fourier_transform_for_square_lattice(k1, k2, unit_cell=0, hopping_1=1, hopping_2=1) + return hamiltonian + +def hamiltonian_of_square_lattice_in_quasi_one_dimension(k, N=10): + 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 + for i in range(N): + h01[i, i] = 1 + hamiltonian = guan.one_dimensional_fourier_transform(k, unit_cell=h00, hopping=h01) + return hamiltonian + +def hamiltonian_of_cubic_lattice(k1, k2, k3): + hamiltonian = guan.three_dimensional_fourier_transform_for_cubic_lattice(k1, k2, k3, unit_cell=0, hopping_1=1, hopping_2=1, hopping_3=1) + return hamiltonian + +def hamiltonian_of_ssh_model(k, v=0.6, w=1): + 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) + return hamiltonian + +def hamiltonian_of_graphene(k1, k2, M=0, t=1, a=1/sqrt(3)): + h0 = np.zeros((2, 2), dtype=complex) # mass term + h1 = np.zeros((2, 2), dtype=complex) # nearest hopping + h0[0, 0] = M + h0[1, 1] = -M + h1[1, 0] = t*(cmath.exp(1j*k2*a)+cmath.exp(1j*sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j/2*k2*a)) + h1[0, 1] = h1[1, 0].conj() + hamiltonian = h0 + h1 + return hamiltonian + +def hamiltonian_of_graphene_with_zigzag_in_quasi_one_dimension(k, N=10, M=0, t=1): + 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 + 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) + return hamiltonian + +def hamiltonian_of_haldane_model(k1, k2, M=2/3, t1=1, t2=1/3, phi=pi/4, a=1/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*sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*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*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a)) + h2[1, 1] = t2*cmath.exp(1j*phi)*(cmath.exp(1j*sqrt(3)*k1*a)+cmath.exp(-1j*sqrt(3)/2*k1*a+1j*3/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j*3/2*k2*a)) + hamiltonian = h0 + h1 + h2 + h2.transpose().conj() + return hamiltonian + +def hamiltonian_of_haldane_model_in_quasi_one_dimension(k, N=10, M=2/3, t1=1, t2=1/3, phi=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() + 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) + return hamiltonian + +def hamiltonian_of_one_QAH_model(k1, k2, t1=1, t2=1, t3=0.5, m=-1): + hamiltonian = np.zeros((2, 2), dtype=complex) + hamiltonian[0, 1] = 2*t1*cos(k1)-1j*2*t1*cos(k2) + hamiltonian[1, 0] = 2*t1*cos(k1)+1j*2*t1*cos(k2) + hamiltonian[0, 0] = m+2*t3*sin(k1)+2*t3*sin(k2)+2*t2*cos(k1+k2) + hamiltonian[1, 1] = -(m+2*t3*sin(k1)+2*t3*sin(k2)+2*t2*cos(k1+k2)) + return hamiltonian + + + + + + + + + + +# Module 5: band_structures_and_wave_functions + +## band structures + +def calculate_eigenvalue(hamiltonian): + if np.array(hamiltonian).shape==(): + eigenvalue = np.real(hamiltonian) + else: + eigenvalue, eigenvector = np.linalg.eigh(hamiltonian) + return eigenvalue + +def calculate_eigenvalue_with_one_parameter(x_array, hamiltonian_function, print_show=0): + 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 + return eigenvalue_array + +def calculate_eigenvalue_with_two_parameters(x_array, y_array, hamiltonian_function, print_show=0, print_show_more=0): + 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 + return eigenvalue_array + +## wave functions + +def calculate_eigenvector(hamiltonian): + eigenvalue, eigenvector = np.linalg.eigh(hamiltonian) + return eigenvector + +## find vector with the same gauge + +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=10001, precision=1e-6): + 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) + return vector_target + +def find_vector_with_fixed_gauge_by_making_one_component_real(vector, precision=0.005, index=None): + 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 + return vector + + + + + + + + + + +# Module 6: Green functions + +def green_function(fermi_energy, hamiltonian, broadening, self_energy=0): + 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) + return green + +def green_function_nn_n(fermi_energy, h00, h01, green_nn_n_minus, broadening, self_energy=0): + 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) + return green_nn_n + +def green_function_in_n(green_in_n_minus, h01, green_nn_n): + green_in_n = np.dot(np.dot(green_in_n_minus, h01), green_nn_n) + return green_in_n + +def green_function_ni_n(green_nn_n, h01, green_ni_n_minus): + h01 = np.array(h01) + green_ni_n = np.dot(np.dot(green_nn_n, h01.transpose().conj()), green_ni_n_minus) + return green_ni_n + +def green_function_ii_n(green_ii_n_minus, green_in_n_minus, h01, green_nn_n, green_ni_n_minus): + 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) + return green_ii_n + +def transfer_matrix(fermi_energy, h00, h01): + 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 + return transfer + +def surface_green_function_of_lead(fermi_energy, h00, h01): + 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))) + return right_lead_surface, left_lead_surface + +def self_energy_of_lead(fermi_energy, h00, h01): + 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 + 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): + 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 + 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): + 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 + return self_energy, gamma + +def green_function_with_leads(fermi_energy, h00, h01, h_LC, h_CR, center_hamiltonian): + 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) + return green, gamma_right, gamma_left + +def electron_correlation_function_green_n_for_local_current(fermi_energy, h00, h01, h_LC, h_CR, center_hamiltonian): + 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())) + return G_n + + + + + + + + + + +# Module 7: density of states + +def total_density_of_states(fermi_energy, hamiltonian, broadening=0.01): + green = guan.green_function(fermi_energy, hamiltonian, broadening) + total_dos = -np.trace(np.imag(green))/pi + return total_dos + +def total_density_of_states_with_fermi_energy_array(fermi_energy_array, hamiltonian, broadening=0.01): + dim = np.array(fermi_energy_array).shape[0] + total_dos_array = np.zeros(dim) + i0 = 0 + for fermi_energy in fermi_energy_array: + total_dos_array[i0] = total_density_of_states(fermi_energy, hamiltonian, broadening) + i0 += 1 + return total_dos_array + +def local_density_of_states_for_square_lattice(fermi_energy, hamiltonian, N1, N2, internal_degree=1, broadening=0.01): + # dim_hamiltonian = N1*N2*internal_degree + 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])/pi + return local_dos + +def local_density_of_states_for_cubic_lattice(fermi_energy, hamiltonian, N1, N2, N3, internal_degree=1, broadening=0.01): + # dim_hamiltonian = N1*N2*N3*internal_degree + 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])/pi + return local_dos + +def local_density_of_states_for_square_lattice_using_dyson_equation(fermi_energy, h00, h01, N2, N1, internal_degree=1, broadening=0.01): + # dim_h00 = N2*internal_degree + 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])/pi + return local_dos + +def local_density_of_states_for_cubic_lattice_using_dyson_equation(fermi_energy, h00, h01, N3, N2, N1, internal_degree=1, broadening=0.01): + # dim_h00 = N2*N3*internal_degree + 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])/pi + return local_dos + +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): + # dim_h00 = N2*internal_degree + 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])/pi + return local_dos + + + + + + + + + + +# Module 8: quantum transport + +## conductance + +def calculate_conductance(fermi_energy, h00, h01, length=100): + right_self_energy, left_self_energy, gamma_right, gamma_left = guan.self_energy_of_lead(fermi_energy, h00, h01) + for ix in range(length): + if ix == 0: + green_nn_n = guan.green_function(fermi_energy, h00, broadening=0, self_energy=left_self_energy) + green_0n_n = copy.deepcopy(green_nn_n) + elif ix != length-1: + green_nn_n = guan.green_function_nn_n(fermi_energy, h00, h01, green_nn_n, broadening=0) + green_0n_n = guan.green_function_in_n(green_0n_n, h01, green_nn_n) + else: + green_nn_n = guan.green_function_nn_n(fermi_energy, h00, h01, green_nn_n, broadening=0, self_energy=right_self_energy) + green_0n_n = guan.green_function_in_n(green_0n_n, h01, green_nn_n) + conductance = np.trace(np.dot(np.dot(np.dot(gamma_left, green_0n_n), gamma_right), green_0n_n.transpose().conj())) + return conductance + +def calculate_conductance_with_fermi_energy_array(fermi_energy_array, h00, h01, length=100): + dim = np.array(fermi_energy_array).shape[0] + conductance_array = np.zeros(dim) + i0 = 0 + for fermi_energy_0 in fermi_energy_array: + conductance_array[i0] = np.real(calculate_conductance(fermi_energy_0, h00, h01, length)) + i0 += 1 + return conductance_array + +def calculate_conductance_with_disorder(fermi_energy, h00, h01, disorder_intensity=2.0, disorder_concentration=1.0, length=100): + right_self_energy, left_self_energy, gamma_right, gamma_left = guan.self_energy_of_lead(fermi_energy, h00, h01) + dim = np.array(h00).shape[0] + for ix 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 ix == 0: + green_nn_n = guan.green_function(fermi_energy, h00+disorder, broadening=0, self_energy=left_self_energy) + green_0n_n = copy.deepcopy(green_nn_n) + elif ix != length-1: + green_nn_n = guan.green_function_nn_n(fermi_energy, h00+disorder, h01, green_nn_n, broadening=0) + green_0n_n = guan.green_function_in_n(green_0n_n, h01, green_nn_n) + else: + green_nn_n = guan.green_function_nn_n(fermi_energy, h00+disorder, h01, green_nn_n, broadening=0, self_energy=right_self_energy) + green_0n_n = guan.green_function_in_n(green_0n_n, h01, green_nn_n) + conductance = np.trace(np.dot(np.dot(np.dot(gamma_left, green_0n_n), gamma_right), green_0n_n.transpose().conj())) + return conductance + +def calculate_conductance_with_disorder_intensity_array(fermi_energy, h00, h01, disorder_intensity_array, disorder_concentration=1.0, length=100, calculation_times=1): + dim = np.array(disorder_intensity_array).shape[0] + conductance_array = np.zeros(dim) + i0 = 0 + for disorder_intensity_0 in disorder_intensity_array: + for times in range(calculation_times): + conductance_array[i0] = conductance_array[i0]+np.real(calculate_conductance_with_disorder(fermi_energy, h00, h01, disorder_intensity=disorder_intensity_0, disorder_concentration=disorder_concentration, length=length)) + i0 += 1 + conductance_array = conductance_array/calculation_times + return conductance_array + +def calculate_conductance_with_disorder_concentration_array(fermi_energy, h00, h01, disorder_concentration_array, disorder_intensity=2.0, length=100, calculation_times=1): + dim = np.array(disorder_concentration_array).shape[0] + conductance_array = np.zeros(dim) + i0 = 0 + for disorder_concentration_0 in disorder_concentration_array: + for times in range(calculation_times): + conductance_array[i0] = conductance_array[i0]+np.real(calculate_conductance_with_disorder(fermi_energy, h00, h01, disorder_intensity=disorder_intensity, disorder_concentration=disorder_concentration_0, length=length)) + i0 += 1 + conductance_array = conductance_array/calculation_times + return conductance_array + +def calculate_conductance_with_scattering_length_array(fermi_energy, h00, h01, length_array, disorder_intensity=2.0, disorder_concentration=1.0, calculation_times=1): + dim = np.array(length_array).shape[0] + conductance_array = np.zeros(dim) + i0 = 0 + for length_0 in length_array: + for times in range(calculation_times): + conductance_array[i0] = conductance_array[i0]+np.real(calculate_conductance_with_disorder(fermi_energy, h00, h01, disorder_intensity=disorder_intensity, disorder_concentration=disorder_concentration, length=length_0)) + i0 += 1 + conductance_array = conductance_array/calculation_times + return conductance_array + + +## scattering matrix + +def if_active_channel(k_of_channel): + if np.abs(np.imag(k_of_channel))<1e-6: + if_active = 1 + else: + if_active = 0 + return if_active + +def get_k_and_velocity_of_channel(fermi_energy, h00, h01): + if np.array(h00).shape==(): + dim = 1 + else: + dim = np.array(h00).shape[0] + transfer = guan.transfer_matrix(fermi_energy, h00, h01) + eigenvalue, eigenvector = np.linalg.eig(transfer) + k_of_channel = np.log(eigenvalue)/1j + ind = np.argsort(np.real(k_of_channel)) + k_of_channel = np.sort(k_of_channel) + temp = np.zeros((2*dim, 2*dim), dtype=complex) + temp2 = np.zeros((2*dim), dtype=complex) + i0 = 0 + for ind0 in ind: + temp[:, i0] = eigenvector[:, ind0] + temp2[i0] = eigenvalue[ind0] + i0 += 1 + eigenvalue = copy.deepcopy(temp2) + temp = temp[0:dim, :] + factor = np.zeros(2*dim, dtype=complex) + for dim0 in range(dim): + factor = factor+np.square(np.abs(temp[dim0, :])) + for dim0 in range(2*dim): + temp[:, dim0] = temp[:, dim0]/np.sqrt(factor[dim0]) + velocity_of_channel = np.zeros((2*dim), dtype=complex) + for dim0 in range(2*dim): + velocity_of_channel[dim0] = eigenvalue[dim0]*np.dot(np.dot(temp[0:dim, :].transpose().conj(), h01),temp[0:dim, :])[dim0, dim0] + velocity_of_channel = -2*np.imag(velocity_of_channel) + eigenvector = copy.deepcopy(temp) + return k_of_channel, velocity_of_channel, eigenvalue, eigenvector + +def get_classified_k_velocity_u_and_f(fermi_energy, h00, h01): + if np.array(h00).shape==(): + dim = 1 + else: + dim = np.array(h00).shape[0] + k_of_channel, velocity_of_channel, eigenvalue, eigenvector = get_k_and_velocity_of_channel(fermi_energy, h00, h01) + ind_right_active = 0; ind_right_evanescent = 0; ind_left_active = 0; ind_left_evanescent = 0 + k_right = np.zeros(dim, dtype=complex); k_left = np.zeros(dim, dtype=complex) + velocity_right = np.zeros(dim, dtype=complex); velocity_left = np.zeros(dim, dtype=complex) + lambda_right = np.zeros(dim, dtype=complex); lambda_left = np.zeros(dim, dtype=complex) + u_right = np.zeros((dim, dim), dtype=complex); u_left = np.zeros((dim, dim), dtype=complex) + for dim0 in range(2*dim): + if_active = if_active_channel(k_of_channel[dim0]) + if if_active_channel(k_of_channel[dim0]) == 1: + direction = np.sign(velocity_of_channel[dim0]) + else: + direction = np.sign(np.imag(k_of_channel[dim0])) + if direction == 1: + if if_active == 1: # right-moving active channel + k_right[ind_right_active] = k_of_channel[dim0] + velocity_right[ind_right_active] = velocity_of_channel[dim0] + lambda_right[ind_right_active] = eigenvalue[dim0] + u_right[:, ind_right_active] = eigenvector[:, dim0] + ind_right_active += 1 + else: # right-moving evanescent channel + k_right[dim-1-ind_right_evanescent] = k_of_channel[dim0] + velocity_right[dim-1-ind_right_evanescent] = velocity_of_channel[dim0] + lambda_right[dim-1-ind_right_evanescent] = eigenvalue[dim0] + u_right[:, dim-1-ind_right_evanescent] = eigenvector[:, dim0] + ind_right_evanescent += 1 + else: + if if_active == 1: # left-moving active channel + k_left[ind_left_active] = k_of_channel[dim0] + velocity_left[ind_left_active] = velocity_of_channel[dim0] + lambda_left[ind_left_active] = eigenvalue[dim0] + u_left[:, ind_left_active] = eigenvector[:, dim0] + ind_left_active += 1 + else: # left-moving evanescent channel + k_left[dim-1-ind_left_evanescent] = k_of_channel[dim0] + velocity_left[dim-1-ind_left_evanescent] = velocity_of_channel[dim0] + lambda_left[dim-1-ind_left_evanescent] = eigenvalue[dim0] + u_left[:, dim-1-ind_left_evanescent] = eigenvector[:, dim0] + ind_left_evanescent += 1 + lambda_matrix_right = np.diag(lambda_right) + lambda_matrix_left = np.diag(lambda_left) + f_right = np.dot(np.dot(u_right, lambda_matrix_right), np.linalg.inv(u_right)) + f_left = np.dot(np.dot(u_left, lambda_matrix_left), np.linalg.inv(u_left)) + return k_right, k_left, velocity_right, velocity_left, f_right, f_left, u_right, u_left, ind_right_active + +def calculate_scattering_matrix(fermi_energy, h00, h01, length=100): + 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 = 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): + 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, 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 = if_active_channel(k_right[dim0])*if_active_channel(k_right[dim1]) + if if_active == 1: + transmission_matrix[dim0, dim1] = np.sqrt(np.abs(velocity_right[dim0]/velocity_right[dim1])) * transmission_matrix[dim0, dim1] + reflection_matrix[dim0, dim1] = np.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!') + return transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active + +def print_or_write_scattering_matrix(fermi_energy, h00, h01, length=100, on_print=1, on_write=0): + if np.array(h00).shape==(): + dim = 1 + else: + dim = np.array(h00).shape[0] + transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active = calculate_scattering_matrix(fermi_energy, h00, h01, length) + if on_print == 1: + print('\nActive channel (left or right) = ', ind_right_active) + print('Evanescent channel (left or right) = ', dim-ind_right_active, '\n') + print('K of right-moving active channels:\n', np.real(k_right[0:ind_right_active])) + print('K of left-moving active channels:\n', np.real(k_left[0:ind_right_active]), '\n') + print('Velocity of right-moving active channels:\n', np.real(velocity_right[0:ind_right_active])) + print('Velocity of left-moving active channels:\n', np.real(velocity_left[0:ind_right_active]), '\n') + print('Transmission matrix:\n', np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active]))) + print('Reflection matrix:\n', np.square(np.abs(reflection_matrix[0:ind_right_active, 0:ind_right_active])), '\n') + print('Total transmission of channels:\n', np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active])), axis=0)) + print('Total reflection of channels:\n',np.sum(np.square(np.abs(reflection_matrix[0:ind_right_active, 0:ind_right_active])), axis=0)) + print('Sum of transmission and reflection of channels:\n', 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)) + print('Total conductance = ', np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active]))), '\n') + if on_write == 1: + with open('a.txt', 'w') as f: + f.write('Active channel (left or right) = ' + str(ind_right_active) + '\n') + f.write('Evanescent channel (left or right) = ' + str(dim - ind_right_active) + '\n\n') + f.write('Channel K Velocity\n') + for ind0 in range(ind_right_active): + f.write(' '+str(ind0 + 1) + ' | '+str(np.real(k_right[ind0]))+' ' + str(np.real(velocity_right[ind0]))+'\n') + f.write('\n') + for ind0 in range(ind_right_active): + f.write(' -' + str(ind0 + 1) + ' | ' + str(np.real(k_left[ind0])) + ' ' + str(np.real(velocity_left[ind0])) + '\n') + f.write('\nScattering matrix:\n ') + for ind0 in range(ind_right_active): + f.write(str(ind0+1)+' ') + f.write('\n') + for ind1 in range(ind_right_active): + f.write(' '+str(ind1+1)+' ') + for ind2 in range(ind_right_active): + f.write('%f' % np.square(np.abs(transmission_matrix[ind1, ind2]))+' ') + f.write('\n') + f.write('\n') + for ind1 in range(ind_right_active): + f.write(' -'+str(ind1+1)+' ') + for ind2 in range(ind_right_active): + f.write('%f' % np.square(np.abs(reflection_matrix[ind1, ind2]))+' ') + f.write('\n') + f.write('\n') + f.write('Total transmission of channels:\n'+str(np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active])), axis=0))+'\n') + f.write('Total conductance = '+str(np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active]))))+'\n') + + + + + + + + + + +# Module 9: topological invariant + +def calculate_chern_number_for_square_lattice(hamiltonian_function, precision=100): + if np.array(hamiltonian_function(0, 0)).shape==(): + dim = 1 + else: + dim = np.array(hamiltonian_function(0, 0)).shape[0] + delta = 2*pi/precision + chern_number = np.zeros(dim, dtype=complex) + for kx in np.arange(-pi, pi, delta): + for ky in np.arange(-pi, 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*pi*1j) + return chern_number + +def calculate_chern_number_for_square_lattice_with_Wilson_loop(hamiltonian_function, precision_of_plaquettes=10, precision_of_Wilson_loop=100): + delta = 2*pi/precision_of_plaquettes + chern_number = 0 + for kx in np.arange(-pi, pi, delta): + for ky in np.arange(-pi, pi, delta): + vector_array = [] + # line_1 + for i0 in range(precision_of_Wilson_loop+1): + 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+1)) + 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+1), 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-1): + H_delta = hamiltonian_function(kx, ky+delta-delta/precision_of_Wilson_loop*(i0+1)) + 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*pi) + return chern_number + +def calculate_chern_number_for_honeycomb_lattice(hamiltonian_function, a=1, precision=300): + 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*sqrt(3)*pi/9/a + L2 = 2*sqrt(3)*pi/9/a + L3 = 2*pi/3/a + delta1 = 2*L1/precision + delta3 = 2*L3/precision + for kx in np.arange(-L1, L1, delta1): + for ky in np.arange(-L3, L3, delta3): + if (-L2<=kx<=L2) or (kx>L2 and -(L1-kx)*tan(pi/3)<=ky<=(L1-kx)*tan(pi/3)) or (kx<-L2 and -(kx-(-L1))*tan(pi/3)<=ky<=(kx-(-L1))*tan(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*pi*1j) + return chern_number + +def calculate_wilson_loop(hamiltonian_function, k_min=-pi, k_max=pi, precision=100): + 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): + 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]) + return wilson_loop_array + + + + + + + + + + +# Module 10: read and write + +def read_one_dimensional_data(filename='a', format='txt'): + f = open(filename+'.'+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) + return x_array, y_array + +def read_two_dimensional_data(filename='a', format='txt'): + f = open(filename+'.'+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) + return x_array, y_array, matrix + +def write_one_dimensional_data(x_array, y_array, filename='a', format='txt'): + x_array = np.array(x_array) + y_array = np.array(y_array) + with open(filename+'.'+format, 'w') as f: + 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 + +def write_two_dimensional_data(x_array, y_array, matrix, filename='a', format='txt'): + x_array = np.array(x_array) + y_array = np.array(y_array) + matrix = np.array(matrix) + with open(filename+'.'+format, 'w') as f: + 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 + + + + + + + + + + +# Module 11: plot figures + +def plot(x_array, y_array, xlabel='x', ylabel='y', title='', show=1, save=0, filename='a', format='jpg', dpi=300, type='', y_min=None, y_max=None, linewidth=None, markersize=None): + import matplotlib.pyplot as plt + fig, ax = plt.subplots() + plt.subplots_adjust(bottom=0.20, left=0.18) + ax.plot(x_array, y_array, type, linewidth=linewidth, markersize=markersize) + ax.grid() + ax.set_title(title, fontsize=20, fontfamily='Times New Roman') + ax.set_xlabel(xlabel, fontsize=20, fontfamily='Times New Roman') + ax.set_ylabel(ylabel, fontsize=20, 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) + ax.tick_params(labelsize=20) + labels = ax.get_xticklabels() + ax.get_yticklabels() + [label.set_fontname('Times New Roman') for label in labels] + if save == 1: + plt.savefig(filename+'.'+format, dpi=dpi) + if show == 1: + plt.show() + plt.close('all') + +def plot_3d_surface(x_array, y_array, matrix, xlabel='x', ylabel='y', zlabel='z', title='', show=1, save=0, filename='a', format='jpg', dpi=300, z_min=None, z_max=None, rcount=100, ccount=100): + 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=20, fontfamily='Times New Roman') + ax.set_xlabel(xlabel, fontsize=20, fontfamily='Times New Roman') + ax.set_ylabel(ylabel, fontsize=20, fontfamily='Times New Roman') + ax.set_zlabel(zlabel, fontsize=20, 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=15) + labels = ax.get_xticklabels() + ax.get_yticklabels() + ax.get_zticklabels() + [label.set_fontname('Times New Roman') for label in labels] + cax = plt.axes([0.80, 0.15, 0.05, 0.75]) + cbar = fig.colorbar(surf, cax=cax) + cbar.ax.tick_params(labelsize=15) + for l in cbar.ax.yaxis.get_ticklabels(): + l.set_family('Times New Roman') + if save == 1: + plt.savefig(filename+'.'+format, dpi=dpi) + if show == 1: + plt.show() + plt.close('all') + +def plot_contour(x_array, y_array, matrix, xlabel='x', ylabel='y', title='', show=1, save=0, filename='a', format='jpg', dpi=300): + import matplotlib.pyplot as plt + fig, ax = plt.subplots() + plt.subplots_adjust(bottom=0.2, right=0.75, left = 0.16) + x_array, y_array = np.meshgrid(x_array, y_array) + contour = ax.contourf(x_array,y_array,matrix,cmap='jet') + ax.set_title(title, fontsize=20, fontfamily='Times New Roman') + ax.set_xlabel(xlabel, fontsize=20, fontfamily='Times New Roman') + ax.set_ylabel(ylabel, fontsize=20, fontfamily='Times New Roman') + ax.tick_params(labelsize=15) + labels = ax.get_xticklabels() + ax.get_yticklabels() + [label.set_fontname('Times New Roman') for label in labels] + cax = plt.axes([0.78, 0.17, 0.08, 0.71]) + cbar = fig.colorbar(contour, cax=cax) + cbar.ax.tick_params(labelsize=15) + for l in cbar.ax.yaxis.get_ticklabels(): + l.set_family('Times New Roman') + if save == 1: + plt.savefig(filename+'.'+format, dpi=dpi) + if show == 1: + plt.show() + plt.close('all') + + + + + + + + + + +# Module 12: others + +## download + +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.iframe['src'] + 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') + +## audio + +def str_to_audio(str='hello world', rate=125, voice=1, read=1, save=0, print_text=0): + import pyttsx3 + if print_text==1: + print(str) + engine = pyttsx3.init() + voices = engine.getProperty('voices') + engine.setProperty('voice', voices[voice].id) + engine.setProperty("rate", rate) + if save==1: + engine.save_to_file(str, 'str.mp3') + engine.runAndWait() + print('MP3 file saved!') + if read==1: + engine.say(str) + engine.runAndWait() + +def txt_to_audio(txt_path, rate=125, voice=1, read=1, save=0, print_text=0): + import pyttsx3 + f = open(txt_path, 'r', encoding ='utf-8') + text = f.read() + if print_text==1: + print(text) + engine = pyttsx3.init() + voices = engine.getProperty('voices') + engine.setProperty('voice', voices[voice].id) + engine.setProperty("rate", rate) + if save==1: + import re + file_name = re.split('[/,\\\]', txt_path)[-1][:-4] + engine.save_to_file(text, file_name+'.mp3') + engine.runAndWait() + print('MP3 file saved!') + if read==1: + engine.say(text) + engine.runAndWait() + +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() + return content + +def pdf_to_audio(pdf_path, rate=125, voice=1, read=1, save=0, print_text=0): + import pyttsx3 + text = guan.pdf_to_text(pdf_path) + text = text.replace('\n', ' ') + if print_text==1: + print(text) + engine = pyttsx3.init() + voices = engine.getProperty('voices') + engine.setProperty('voice', voices[voice].id) + engine.setProperty("rate", rate) + if save==1: + import re + file_name = re.split('[/,\\\]', pdf_path)[-1][:-4] + engine.save_to_file(text, file_name+'.mp3') + engine.runAndWait() + print('MP3 file saved!') + if read==1: + engine.say(text) + engine.runAndWait() \ 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 deleted file mode 100644 index 480bec6..0000000 --- a/PyPI/src/guan/band_structures_and_wave_functions.py +++ /dev/null @@ -1,125 +0,0 @@ -# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about. The primary location of this package is on website https://py.guanjihuan.com. - -# band_structures_and_wave_functions - -## band structures - -import numpy as np -import cmath - -def calculate_eigenvalue(hamiltonian): - if np.array(hamiltonian).shape==(): - eigenvalue = np.real(hamiltonian) - else: - eigenvalue, eigenvector = np.linalg.eigh(hamiltonian) - return eigenvalue - -def calculate_eigenvalue_with_one_parameter(x_array, hamiltonian_function, print_show=0): - 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 - return eigenvalue_array - -def calculate_eigenvalue_with_two_parameters(x_array, y_array, hamiltonian_function, print_show=0, print_show_more=0): - 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 - return eigenvalue_array - -## wave functions - -def calculate_eigenvector(hamiltonian): - eigenvalue, eigenvector = np.linalg.eigh(hamiltonian) - return eigenvector - -## find vector with the same gauge - -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=10001, precision=1e-6): - 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) - return vector_target - -def find_vector_with_fixed_gauge_by_making_one_component_real(vector, precision=0.005, index=None): - 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 - return vector \ No newline at end of file diff --git a/PyPI/src/guan/basic_functions.py b/PyPI/src/guan/basic_functions.py deleted file mode 100644 index 84bb571..0000000 --- a/PyPI/src/guan/basic_functions.py +++ /dev/null @@ -1,74 +0,0 @@ -# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about. The primary location of this package is on website https://py.guanjihuan.com. - -# basic functions - -import numpy as np - -## test - -def test(): - print('\nSuccess in the installation of Guan package!\n') - -## Pauli matrices - -def sigma_0(): - return np.eye(2) - -def sigma_x(): - return np.array([[0, 1],[1, 0]]) - -def sigma_y(): - return np.array([[0, -1j],[1j, 0]]) - -def sigma_z(): - return np.array([[1, 0],[0, -1]]) - -## Kronecker product of Pauli matrices - -def sigma_00(): - return np.kron(sigma_0(), sigma_0()) - -def sigma_0x(): - return np.kron(sigma_0(), sigma_x()) - -def sigma_0y(): - return np.kron(sigma_0(), sigma_y()) - -def sigma_0z(): - return np.kron(sigma_0(), sigma_z()) - -def sigma_x0(): - return np.kron(sigma_x(), sigma_0()) - -def sigma_xx(): - return np.kron(sigma_x(), sigma_x()) - -def sigma_xy(): - return np.kron(sigma_x(), sigma_y()) - -def sigma_xz(): - return np.kron(sigma_x(), sigma_z()) - -def sigma_y0(): - return np.kron(sigma_y(), sigma_0()) - -def sigma_yx(): - return np.kron(sigma_y(), sigma_x()) - -def sigma_yy(): - return np.kron(sigma_y(), sigma_y()) - -def sigma_yz(): - return np.kron(sigma_y(), sigma_z()) - -def sigma_z0(): - return np.kron(sigma_z(), sigma_0()) - -def sigma_zx(): - return np.kron(sigma_z(), sigma_x()) - -def sigma_zy(): - return np.kron(sigma_z(), sigma_y()) - -def sigma_zz(): - return np.kron(sigma_z(), sigma_z()) \ No newline at end of file diff --git a/PyPI/src/guan/density_of_states.py b/PyPI/src/guan/density_of_states.py deleted file mode 100644 index 7065bf2..0000000 --- a/PyPI/src/guan/density_of_states.py +++ /dev/null @@ -1,139 +0,0 @@ -# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about. The primary location of this package is on website https://py.guanjihuan.com. - -# density of states - -import numpy as np -from math import * -import guan - -def total_density_of_states(fermi_energy, hamiltonian, broadening=0.01): - green = guan.green_function(fermi_energy, hamiltonian, broadening) - total_dos = -np.trace(np.imag(green))/pi - return total_dos - -def total_density_of_states_with_fermi_energy_array(fermi_energy_array, hamiltonian, broadening=0.01): - dim = np.array(fermi_energy_array).shape[0] - total_dos_array = np.zeros(dim) - i0 = 0 - for fermi_energy in fermi_energy_array: - total_dos_array[i0] = total_density_of_states(fermi_energy, hamiltonian, broadening) - i0 += 1 - return total_dos_array - -def local_density_of_states_for_square_lattice(fermi_energy, hamiltonian, N1, N2, internal_degree=1, broadening=0.01): - # dim_hamiltonian = N1*N2*internal_degree - 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])/pi - return local_dos - -def local_density_of_states_for_cubic_lattice(fermi_energy, hamiltonian, N1, N2, N3, internal_degree=1, broadening=0.01): - # dim_hamiltonian = N1*N2*N3*internal_degree - 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])/pi - return local_dos - -def local_density_of_states_for_square_lattice_using_dyson_equation(fermi_energy, h00, h01, N2, N1, internal_degree=1, broadening=0.01): - # dim_h00 = N2*internal_degree - 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])/pi - return local_dos - -def local_density_of_states_for_cubic_lattice_using_dyson_equation(fermi_energy, h00, h01, N3, N2, N1, internal_degree=1, broadening=0.01): - # dim_h00 = N2*N3*internal_degree - 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])/pi - return local_dos - -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): - # dim_h00 = N2*internal_degree - 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])/pi - return local_dos \ No newline at end of file diff --git a/PyPI/src/guan/others.py b/PyPI/src/guan/others.py deleted file mode 100644 index 7ae9620..0000000 --- a/PyPI/src/guan/others.py +++ /dev/null @@ -1,130 +0,0 @@ -# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about. The primary location of this package is on website https://py.guanjihuan.com. - -# others - -import guan - -## download - -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.iframe['src'] - 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') - -## audio - -def str_to_audio(str='hello world', rate=125, voice=1, read=1, save=0, print_text=0): - import pyttsx3 - if print_text==1: - print(str) - engine = pyttsx3.init() - voices = engine.getProperty('voices') - engine.setProperty('voice', voices[voice].id) - engine.setProperty("rate", rate) - if save==1: - engine.save_to_file(str, 'str.mp3') - engine.runAndWait() - print('MP3 file saved!') - if read==1: - engine.say(str) - engine.runAndWait() - -def txt_to_audio(txt_path, rate=125, voice=1, read=1, save=0, print_text=0): - import pyttsx3 - f = open(txt_path, 'r', encoding ='utf-8') - text = f.read() - if print_text==1: - print(text) - engine = pyttsx3.init() - voices = engine.getProperty('voices') - engine.setProperty('voice', voices[voice].id) - engine.setProperty("rate", rate) - if save==1: - import re - file_name = re.split('[/,\\\]', txt_path)[-1][:-4] - engine.save_to_file(text, file_name+'.mp3') - engine.runAndWait() - print('MP3 file saved!') - if read==1: - engine.say(text) - engine.runAndWait() - -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() - return content - -def pdf_to_audio(pdf_path, rate=125, voice=1, read=1, save=0, print_text=0): - import pyttsx3 - text = guan.pdf_to_text(pdf_path) - text = text.replace('\n', ' ') - if print_text==1: - print(text) - engine = pyttsx3.init() - voices = engine.getProperty('voices') - engine.setProperty('voice', voices[voice].id) - engine.setProperty("rate", rate) - if save==1: - import re - file_name = re.split('[/,\\\]', pdf_path)[-1][:-4] - engine.save_to_file(text, file_name+'.mp3') - engine.runAndWait() - print('MP3 file saved!') - if read==1: - engine.say(text) - engine.runAndWait() \ No newline at end of file diff --git a/PyPI/src/guan/plot_figures.py b/PyPI/src/guan/plot_figures.py deleted file mode 100644 index b35b2d5..0000000 --- a/PyPI/src/guan/plot_figures.py +++ /dev/null @@ -1,91 +0,0 @@ -# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about. The primary location of this package is on website https://py.guanjihuan.com. - -# plot figures - -import numpy as np - -def plot(x_array, y_array, xlabel='x', ylabel='y', title='', show=1, save=0, filename='a', format='jpg', dpi=300, type='', y_min=None, y_max=None, linewidth=None, markersize=None): - import matplotlib.pyplot as plt - fig, ax = plt.subplots() - plt.subplots_adjust(bottom=0.20, left=0.18) - ax.plot(x_array, y_array, type, linewidth=linewidth, markersize=markersize) - ax.grid() - ax.set_title(title, fontsize=20, fontfamily='Times New Roman') - ax.set_xlabel(xlabel, fontsize=20, fontfamily='Times New Roman') - ax.set_ylabel(ylabel, fontsize=20, 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) - ax.tick_params(labelsize=20) - labels = ax.get_xticklabels() + ax.get_yticklabels() - [label.set_fontname('Times New Roman') for label in labels] - if save == 1: - plt.savefig(filename+'.'+format, dpi=dpi) - if show == 1: - plt.show() - plt.close('all') - -def plot_3d_surface(x_array, y_array, matrix, xlabel='x', ylabel='y', zlabel='z', title='', show=1, save=0, filename='a', format='jpg', dpi=300, z_min=None, z_max=None, rcount=100, ccount=100): - 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=20, fontfamily='Times New Roman') - ax.set_xlabel(xlabel, fontsize=20, fontfamily='Times New Roman') - ax.set_ylabel(ylabel, fontsize=20, fontfamily='Times New Roman') - ax.set_zlabel(zlabel, fontsize=20, 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=15) - labels = ax.get_xticklabels() + ax.get_yticklabels() + ax.get_zticklabels() - [label.set_fontname('Times New Roman') for label in labels] - cax = plt.axes([0.80, 0.15, 0.05, 0.75]) - cbar = fig.colorbar(surf, cax=cax) - cbar.ax.tick_params(labelsize=15) - for l in cbar.ax.yaxis.get_ticklabels(): - l.set_family('Times New Roman') - if save == 1: - plt.savefig(filename+'.'+format, dpi=dpi) - if show == 1: - plt.show() - plt.close('all') - -def plot_contour(x_array, y_array, matrix, xlabel='x', ylabel='y', title='', show=1, save=0, filename='a', format='jpg', dpi=300): - import matplotlib.pyplot as plt - fig, ax = plt.subplots() - plt.subplots_adjust(bottom=0.2, right=0.75, left = 0.16) - x_array, y_array = np.meshgrid(x_array, y_array) - contour = ax.contourf(x_array,y_array,matrix,cmap='jet') - ax.set_title(title, fontsize=20, fontfamily='Times New Roman') - ax.set_xlabel(xlabel, fontsize=20, fontfamily='Times New Roman') - ax.set_ylabel(ylabel, fontsize=20, fontfamily='Times New Roman') - ax.tick_params(labelsize=15) - labels = ax.get_xticklabels() + ax.get_yticklabels() - [label.set_fontname('Times New Roman') for label in labels] - cax = plt.axes([0.78, 0.17, 0.08, 0.71]) - cbar = fig.colorbar(contour, cax=cax) - cbar.ax.tick_params(labelsize=15) - for l in cbar.ax.yaxis.get_ticklabels(): - l.set_family('Times New Roman') - if save == 1: - plt.savefig(filename+'.'+format, dpi=dpi) - if show == 1: - plt.show() - plt.close('all') \ No newline at end of file diff --git a/PyPI/src/guan/quantum_transport.py b/PyPI/src/guan/quantum_transport.py deleted file mode 100644 index 50e5eab..0000000 --- a/PyPI/src/guan/quantum_transport.py +++ /dev/null @@ -1,264 +0,0 @@ -# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about. The primary location of this package is on website https://py.guanjihuan.com. - -# quantum transport - -## conductance - -import numpy as np -import copy -import guan - -def calculate_conductance(fermi_energy, h00, h01, length=100): - right_self_energy, left_self_energy, gamma_right, gamma_left = guan.self_energy_of_lead(fermi_energy, h00, h01) - for ix in range(length): - if ix == 0: - green_nn_n = guan.green_function(fermi_energy, h00, broadening=0, self_energy=left_self_energy) - green_0n_n = copy.deepcopy(green_nn_n) - elif ix != length-1: - green_nn_n = guan.green_function_nn_n(fermi_energy, h00, h01, green_nn_n, broadening=0) - green_0n_n = guan.green_function_in_n(green_0n_n, h01, green_nn_n) - else: - green_nn_n = guan.green_function_nn_n(fermi_energy, h00, h01, green_nn_n, broadening=0, self_energy=right_self_energy) - green_0n_n = guan.green_function_in_n(green_0n_n, h01, green_nn_n) - conductance = np.trace(np.dot(np.dot(np.dot(gamma_left, green_0n_n), gamma_right), green_0n_n.transpose().conj())) - return conductance - -def calculate_conductance_with_fermi_energy_array(fermi_energy_array, h00, h01, length=100): - dim = np.array(fermi_energy_array).shape[0] - conductance_array = np.zeros(dim) - i0 = 0 - for fermi_energy_0 in fermi_energy_array: - conductance_array[i0] = np.real(calculate_conductance(fermi_energy_0, h00, h01, length)) - i0 += 1 - return conductance_array - -def calculate_conductance_with_disorder(fermi_energy, h00, h01, disorder_intensity=2.0, disorder_concentration=1.0, length=100): - right_self_energy, left_self_energy, gamma_right, gamma_left = guan.self_energy_of_lead(fermi_energy, h00, h01) - dim = np.array(h00).shape[0] - for ix 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 ix == 0: - green_nn_n = guan.green_function(fermi_energy, h00+disorder, broadening=0, self_energy=left_self_energy) - green_0n_n = copy.deepcopy(green_nn_n) - elif ix != length-1: - green_nn_n = guan.green_function_nn_n(fermi_energy, h00+disorder, h01, green_nn_n, broadening=0) - green_0n_n = guan.green_function_in_n(green_0n_n, h01, green_nn_n) - else: - green_nn_n = guan.green_function_nn_n(fermi_energy, h00+disorder, h01, green_nn_n, broadening=0, self_energy=right_self_energy) - green_0n_n = guan.green_function_in_n(green_0n_n, h01, green_nn_n) - conductance = np.trace(np.dot(np.dot(np.dot(gamma_left, green_0n_n), gamma_right), green_0n_n.transpose().conj())) - return conductance - -def calculate_conductance_with_disorder_intensity_array(fermi_energy, h00, h01, disorder_intensity_array, disorder_concentration=1.0, length=100, calculation_times=1): - dim = np.array(disorder_intensity_array).shape[0] - conductance_array = np.zeros(dim) - i0 = 0 - for disorder_intensity_0 in disorder_intensity_array: - for times in range(calculation_times): - conductance_array[i0] = conductance_array[i0]+np.real(calculate_conductance_with_disorder(fermi_energy, h00, h01, disorder_intensity=disorder_intensity_0, disorder_concentration=disorder_concentration, length=length)) - i0 += 1 - conductance_array = conductance_array/calculation_times - return conductance_array - -def calculate_conductance_with_disorder_concentration_array(fermi_energy, h00, h01, disorder_concentration_array, disorder_intensity=2.0, length=100, calculation_times=1): - dim = np.array(disorder_concentration_array).shape[0] - conductance_array = np.zeros(dim) - i0 = 0 - for disorder_concentration_0 in disorder_concentration_array: - for times in range(calculation_times): - conductance_array[i0] = conductance_array[i0]+np.real(calculate_conductance_with_disorder(fermi_energy, h00, h01, disorder_intensity=disorder_intensity, disorder_concentration=disorder_concentration_0, length=length)) - i0 += 1 - conductance_array = conductance_array/calculation_times - return conductance_array - -def calculate_conductance_with_scattering_length_array(fermi_energy, h00, h01, length_array, disorder_intensity=2.0, disorder_concentration=1.0, calculation_times=1): - dim = np.array(length_array).shape[0] - conductance_array = np.zeros(dim) - i0 = 0 - for length_0 in length_array: - for times in range(calculation_times): - conductance_array[i0] = conductance_array[i0]+np.real(calculate_conductance_with_disorder(fermi_energy, h00, h01, disorder_intensity=disorder_intensity, disorder_concentration=disorder_concentration, length=length_0)) - i0 += 1 - conductance_array = conductance_array/calculation_times - return conductance_array - - -## scattering matrix - -def if_active_channel(k_of_channel): - if np.abs(np.imag(k_of_channel))<1e-6: - if_active = 1 - else: - if_active = 0 - return if_active - -def get_k_and_velocity_of_channel(fermi_energy, h00, h01): - if np.array(h00).shape==(): - dim = 1 - else: - dim = np.array(h00).shape[0] - transfer = guan.transfer_matrix(fermi_energy, h00, h01) - eigenvalue, eigenvector = np.linalg.eig(transfer) - k_of_channel = np.log(eigenvalue)/1j - ind = np.argsort(np.real(k_of_channel)) - k_of_channel = np.sort(k_of_channel) - temp = np.zeros((2*dim, 2*dim), dtype=complex) - temp2 = np.zeros((2*dim), dtype=complex) - i0 = 0 - for ind0 in ind: - temp[:, i0] = eigenvector[:, ind0] - temp2[i0] = eigenvalue[ind0] - i0 += 1 - eigenvalue = copy.deepcopy(temp2) - temp = temp[0:dim, :] - factor = np.zeros(2*dim, dtype=complex) - for dim0 in range(dim): - factor = factor+np.square(np.abs(temp[dim0, :])) - for dim0 in range(2*dim): - temp[:, dim0] = temp[:, dim0]/np.sqrt(factor[dim0]) - velocity_of_channel = np.zeros((2*dim), dtype=complex) - for dim0 in range(2*dim): - velocity_of_channel[dim0] = eigenvalue[dim0]*np.dot(np.dot(temp[0:dim, :].transpose().conj(), h01),temp[0:dim, :])[dim0, dim0] - velocity_of_channel = -2*np.imag(velocity_of_channel) - eigenvector = copy.deepcopy(temp) - return k_of_channel, velocity_of_channel, eigenvalue, eigenvector - -def get_classified_k_velocity_u_and_f(fermi_energy, h00, h01): - if np.array(h00).shape==(): - dim = 1 - else: - dim = np.array(h00).shape[0] - k_of_channel, velocity_of_channel, eigenvalue, eigenvector = get_k_and_velocity_of_channel(fermi_energy, h00, h01) - ind_right_active = 0; ind_right_evanescent = 0; ind_left_active = 0; ind_left_evanescent = 0 - k_right = np.zeros(dim, dtype=complex); k_left = np.zeros(dim, dtype=complex) - velocity_right = np.zeros(dim, dtype=complex); velocity_left = np.zeros(dim, dtype=complex) - lambda_right = np.zeros(dim, dtype=complex); lambda_left = np.zeros(dim, dtype=complex) - u_right = np.zeros((dim, dim), dtype=complex); u_left = np.zeros((dim, dim), dtype=complex) - for dim0 in range(2*dim): - if_active = if_active_channel(k_of_channel[dim0]) - if if_active_channel(k_of_channel[dim0]) == 1: - direction = np.sign(velocity_of_channel[dim0]) - else: - direction = np.sign(np.imag(k_of_channel[dim0])) - if direction == 1: - if if_active == 1: # right-moving active channel - k_right[ind_right_active] = k_of_channel[dim0] - velocity_right[ind_right_active] = velocity_of_channel[dim0] - lambda_right[ind_right_active] = eigenvalue[dim0] - u_right[:, ind_right_active] = eigenvector[:, dim0] - ind_right_active += 1 - else: # right-moving evanescent channel - k_right[dim-1-ind_right_evanescent] = k_of_channel[dim0] - velocity_right[dim-1-ind_right_evanescent] = velocity_of_channel[dim0] - lambda_right[dim-1-ind_right_evanescent] = eigenvalue[dim0] - u_right[:, dim-1-ind_right_evanescent] = eigenvector[:, dim0] - ind_right_evanescent += 1 - else: - if if_active == 1: # left-moving active channel - k_left[ind_left_active] = k_of_channel[dim0] - velocity_left[ind_left_active] = velocity_of_channel[dim0] - lambda_left[ind_left_active] = eigenvalue[dim0] - u_left[:, ind_left_active] = eigenvector[:, dim0] - ind_left_active += 1 - else: # left-moving evanescent channel - k_left[dim-1-ind_left_evanescent] = k_of_channel[dim0] - velocity_left[dim-1-ind_left_evanescent] = velocity_of_channel[dim0] - lambda_left[dim-1-ind_left_evanescent] = eigenvalue[dim0] - u_left[:, dim-1-ind_left_evanescent] = eigenvector[:, dim0] - ind_left_evanescent += 1 - lambda_matrix_right = np.diag(lambda_right) - lambda_matrix_left = np.diag(lambda_left) - f_right = np.dot(np.dot(u_right, lambda_matrix_right), np.linalg.inv(u_right)) - f_left = np.dot(np.dot(u_left, lambda_matrix_left), np.linalg.inv(u_left)) - return k_right, k_left, velocity_right, velocity_left, f_right, f_left, u_right, u_left, ind_right_active - -def calculate_scattering_matrix(fermi_energy, h00, h01, length=100): - 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 = 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): - 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, 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 = if_active_channel(k_right[dim0])*if_active_channel(k_right[dim1]) - if if_active == 1: - transmission_matrix[dim0, dim1] = np.sqrt(np.abs(velocity_right[dim0]/velocity_right[dim1])) * transmission_matrix[dim0, dim1] - reflection_matrix[dim0, dim1] = np.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!') - return transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active - -def print_or_write_scattering_matrix(fermi_energy, h00, h01, length=100, on_print=1, on_write=0): - if np.array(h00).shape==(): - dim = 1 - else: - dim = np.array(h00).shape[0] - transmission_matrix, reflection_matrix, k_right, k_left, velocity_right, velocity_left, ind_right_active = calculate_scattering_matrix(fermi_energy, h00, h01, length) - if on_print == 1: - print('\nActive channel (left or right) = ', ind_right_active) - print('Evanescent channel (left or right) = ', dim-ind_right_active, '\n') - print('K of right-moving active channels:\n', np.real(k_right[0:ind_right_active])) - print('K of left-moving active channels:\n', np.real(k_left[0:ind_right_active]), '\n') - print('Velocity of right-moving active channels:\n', np.real(velocity_right[0:ind_right_active])) - print('Velocity of left-moving active channels:\n', np.real(velocity_left[0:ind_right_active]), '\n') - print('Transmission matrix:\n', np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active]))) - print('Reflection matrix:\n', np.square(np.abs(reflection_matrix[0:ind_right_active, 0:ind_right_active])), '\n') - print('Total transmission of channels:\n', np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active])), axis=0)) - print('Total reflection of channels:\n',np.sum(np.square(np.abs(reflection_matrix[0:ind_right_active, 0:ind_right_active])), axis=0)) - print('Sum of transmission and reflection of channels:\n', 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)) - print('Total conductance = ', np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active]))), '\n') - if on_write == 1: - with open('a.txt', 'w') as f: - f.write('Active channel (left or right) = ' + str(ind_right_active) + '\n') - f.write('Evanescent channel (left or right) = ' + str(dim - ind_right_active) + '\n\n') - f.write('Channel K Velocity\n') - for ind0 in range(ind_right_active): - f.write(' '+str(ind0 + 1) + ' | '+str(np.real(k_right[ind0]))+' ' + str(np.real(velocity_right[ind0]))+'\n') - f.write('\n') - for ind0 in range(ind_right_active): - f.write(' -' + str(ind0 + 1) + ' | ' + str(np.real(k_left[ind0])) + ' ' + str(np.real(velocity_left[ind0])) + '\n') - f.write('\nScattering matrix:\n ') - for ind0 in range(ind_right_active): - f.write(str(ind0+1)+' ') - f.write('\n') - for ind1 in range(ind_right_active): - f.write(' '+str(ind1+1)+' ') - for ind2 in range(ind_right_active): - f.write('%f' % np.square(np.abs(transmission_matrix[ind1, ind2]))+' ') - f.write('\n') - f.write('\n') - for ind1 in range(ind_right_active): - f.write(' -'+str(ind1+1)+' ') - for ind2 in range(ind_right_active): - f.write('%f' % np.square(np.abs(reflection_matrix[ind1, ind2]))+' ') - f.write('\n') - f.write('\n') - f.write('Total transmission of channels:\n'+str(np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active])), axis=0))+'\n') - f.write('Total conductance = '+str(np.sum(np.square(np.abs(transmission_matrix[0:ind_right_active, 0:ind_right_active]))))+'\n') \ No newline at end of file diff --git a/PyPI/src/guan/read_and_write.py b/PyPI/src/guan/read_and_write.py deleted file mode 100644 index 4aa0886..0000000 --- a/PyPI/src/guan/read_and_write.py +++ /dev/null @@ -1,87 +0,0 @@ -# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about. The primary location of this package is on website https://py.guanjihuan.com. - -# read and write - -import numpy as np - -def read_one_dimensional_data(filename='a', format='txt'): - f = open(filename+'.'+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) - return x_array, y_array - -def read_two_dimensional_data(filename='a', format='txt'): - f = open(filename+'.'+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) - return x_array, y_array, matrix - -def write_one_dimensional_data(x_array, y_array, filename='a', format='txt'): - x_array = np.array(x_array) - y_array = np.array(y_array) - with open(filename+'.'+format, 'w') as f: - 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 - -def write_two_dimensional_data(x_array, y_array, matrix, filename='a', format='txt'): - x_array = np.array(x_array) - y_array = np.array(y_array) - matrix = np.array(matrix) - with open(filename+'.'+format, 'w') as f: - 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 \ No newline at end of file diff --git a/PyPI/src/guan/topological_invariant.py b/PyPI/src/guan/topological_invariant.py deleted file mode 100644 index 221f140..0000000 --- a/PyPI/src/guan/topological_invariant.py +++ /dev/null @@ -1,131 +0,0 @@ -# Guan is an open-source python package developed and maintained by https://www.guanjihuan.com/about. The primary location of this package is on website https://py.guanjihuan.com. - -# topological invariant - -import numpy as np -import cmath -from math import * -import guan - -def calculate_chern_number_for_square_lattice(hamiltonian_function, precision=100): - if np.array(hamiltonian_function(0, 0)).shape==(): - dim = 1 - else: - dim = np.array(hamiltonian_function(0, 0)).shape[0] - delta = 2*pi/precision - chern_number = np.zeros(dim, dtype=complex) - for kx in np.arange(-pi, pi, delta): - for ky in np.arange(-pi, 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*pi*1j) - return chern_number - -def calculate_chern_number_for_square_lattice_with_Wilson_loop(hamiltonian_function, precision_of_plaquettes=10, precision_of_Wilson_loop=100): - delta = 2*pi/precision_of_plaquettes - chern_number = 0 - for kx in np.arange(-pi, pi, delta): - for ky in np.arange(-pi, pi, delta): - vector_array = [] - # line_1 - for i0 in range(precision_of_Wilson_loop+1): - 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+1)) - 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+1), 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-1): - H_delta = hamiltonian_function(kx, ky+delta-delta/precision_of_Wilson_loop*(i0+1)) - 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*pi) - return chern_number - -def calculate_chern_number_for_honeycomb_lattice(hamiltonian_function, a=1, precision=300): - 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*sqrt(3)*pi/9/a - L2 = 2*sqrt(3)*pi/9/a - L3 = 2*pi/3/a - delta1 = 2*L1/precision - delta3 = 2*L3/precision - for kx in np.arange(-L1, L1, delta1): - for ky in np.arange(-L3, L3, delta3): - if (-L2<=kx<=L2) or (kx>L2 and -(L1-kx)*tan(pi/3)<=ky<=(L1-kx)*tan(pi/3)) or (kx<-L2 and -(kx-(-L1))*tan(pi/3)<=ky<=(kx-(-L1))*tan(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*pi*1j) - return chern_number - -def calculate_wilson_loop(hamiltonian_function, k_min=-pi, k_max=pi, precision=100): - 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): - 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]) - return wilson_loop_array \ No newline at end of file diff --git a/Tutorial/01_test_and_Pauli_matrix.py b/Tutorial/01_test_and_Pauli_matrix.py new file mode 100644 index 0000000..8e1edb0 --- /dev/null +++ b/Tutorial/01_test_and_Pauli_matrix.py @@ -0,0 +1,47 @@ +import guan + +# test +guan.test() + +# Pauli matrix +sigma_0 = guan.sigma_0() +sigma_x = guan.sigma_x() +sigma_y = guan.sigma_y() +sigma_z = guan.sigma_z() +sigma_00 = guan.sigma_00() +sigma_0x = guan.sigma_0x() +sigma_0y = guan.sigma_0y() +sigma_0z = guan.sigma_0z() +sigma_x0 = guan.sigma_x0() +sigma_xx = guan.sigma_xx() +sigma_xy = guan.sigma_xy() +sigma_xz = guan.sigma_xz() +sigma_y0 = guan.sigma_y0() +sigma_yx = guan.sigma_yx() +sigma_yy = guan.sigma_yy() +sigma_yz = guan.sigma_yz() +sigma_z0 = guan.sigma_z0() +sigma_zx = guan.sigma_zx() +sigma_zy = guan.sigma_zy() +sigma_zz = guan.sigma_zz() +print('Pauli matrix\n') +print('sigma_0\n', sigma_0, '\n') +print('sigma_x\n', sigma_x, '\n') +print('sigma_y\n', sigma_y, '\n') +print('sigma_z\n', sigma_z, '\n') +print('sigma_00\n', sigma_00, '\n') +print('sigma_0x\n', sigma_0x, '\n') +print('sigma_0y\n', sigma_0y, '\n') +print('sigma_0z\n', sigma_0z, '\n') +print('sigma_x0\n', sigma_x0, '\n') +print('sigma_xx\n', sigma_xx, '\n') +print('sigma_xy\n', sigma_xy, '\n') +print('sigma_xz\n', sigma_xz, '\n') +print('sigma_y0\n', sigma_y0, '\n') +print('sigma_yx\n', sigma_yx, '\n') +print('sigma_yy\n', sigma_yy, '\n') +print('sigma_yz\n', sigma_yz, '\n') +print('sigma_z0\n', sigma_z0, '\n') +print('sigma_zx\n', sigma_zx, '\n') +print('sigma_zy\n', sigma_zy, '\n') +print('sigma_zz\n', sigma_zz, '\n') \ No newline at end of file diff --git a/Tutorial/02_Fourier_transform_and_calculate_band_structures.py.py b/Tutorial/02_Fourier_transform_and_calculate_band_structures.py.py new file mode 100644 index 0000000..64e46da --- /dev/null +++ b/Tutorial/02_Fourier_transform_and_calculate_band_structures.py.py @@ -0,0 +1,8 @@ +import guan +import numpy as np + +# Fourier transform / calculate band structures / plot figures +k_array = np.linspace(-np.pi, np.pi, 100) +hamiltonian_function = guan.one_dimensional_fourier_transform_with_k(unit_cell=0, hopping=1) # one dimensional chain +eigenvalue_array = guan.calculate_eigenvalue_with_one_parameter(k_array, hamiltonian_function) +guan.plot(k_array, eigenvalue_array, xlabel='k', ylabel='E', type='-k') \ No newline at end of file diff --git a/Tutorial/03_Hamiltonian_of_finite_size_systems.py b/Tutorial/03_Hamiltonian_of_finite_size_systems.py new file mode 100644 index 0000000..a3a8afb --- /dev/null +++ b/Tutorial/03_Hamiltonian_of_finite_size_systems.py @@ -0,0 +1,6 @@ +import guan + +# Hamiltonian of finite size +print('\n', guan.hamiltonian_of_finite_size_system_along_one_direction(3), '\n') +print(guan.hamiltonian_of_finite_size_system_along_two_directions_for_square_lattice(2, 2), '\n') +print(guan.hamiltonian_of_finite_size_system_along_three_directions_for_cubic_lattice(2, 2, 2), '\n') \ No newline at end of file diff --git a/Tutorial/04_some_models_in_the_reciprocal_space.py b/Tutorial/04_some_models_in_the_reciprocal_space.py new file mode 100644 index 0000000..554d5b5 --- /dev/null +++ b/Tutorial/04_some_models_in_the_reciprocal_space.py @@ -0,0 +1,9 @@ +import guan +import numpy as np + +# Hamiltonian of models in the reciprocal space / calculate band structures / plot figures +k_array = np.linspace(-np.pi, np.pi, 100) +eigenvalue_array = guan.calculate_eigenvalue_with_one_parameter(k_array, guan.hamiltonian_of_square_lattice_in_quasi_one_dimension) +guan.plot(k_array, eigenvalue_array, xlabel='k', ylabel='E', type='-k') +eigenvalue_array = guan.calculate_eigenvalue_with_one_parameter(k_array, guan.hamiltonian_of_graphene_with_zigzag_in_quasi_one_dimension) +guan.plot(k_array, eigenvalue_array, xlabel='k', ylabel='E', type='-k') \ No newline at end of file diff --git a/Tutorial/calculate_density_of_states.py b/Tutorial/05_calculate_density_of_states.py similarity index 68% rename from Tutorial/calculate_density_of_states.py rename to Tutorial/05_calculate_density_of_states.py index d7160bc..63fef57 100644 --- a/Tutorial/calculate_density_of_states.py +++ b/Tutorial/05_calculate_density_of_states.py @@ -2,7 +2,7 @@ import guan import numpy as np # calculate density of states -hamiltonian = guan.finite_size_along_two_directions_for_square_lattice(2,2) +hamiltonian = guan.hamiltonian_of_finite_size_system_along_two_directions_for_square_lattice(2,2) fermi_energy_array = np.linspace(-4, 4, 400) total_dos_array = guan.total_density_of_states_with_fermi_energy_array(fermi_energy_array, hamiltonian, broadening=0.1) guan.plot(fermi_energy_array, total_dos_array, xlabel='E', ylabel='Total DOS', type='-o') @@ -10,10 +10,10 @@ guan.plot(fermi_energy_array, total_dos_array, xlabel='E', ylabel='Total DOS', t fermi_energy = 0 N1 = 3 N2 = 4 -hamiltonian = guan.finite_size_along_two_directions_for_square_lattice(N1,N2) +hamiltonian = guan.hamiltonian_of_finite_size_system_along_two_directions_for_square_lattice(N1,N2) LDOS = guan.local_density_of_states_for_square_lattice(fermi_energy, hamiltonian, N1=N1, N2=N2) print('square lattice:\n', LDOS, '\n') -h00 = guan.finite_size_along_one_direction(N2) +h00 = guan.hamiltonian_of_finite_size_system_along_one_direction(N2) h01 = np.identity(N2) LDOS = guan.local_density_of_states_for_square_lattice_using_dyson_equation(fermi_energy, h00=h00, h01=h01, N2=N2, N1=N1) print(LDOS, '\n\n') @@ -22,10 +22,10 @@ print(LDOS, '\n\n') N1 = 3 N2 = 4 N3 = 5 -hamiltonian = guan.finite_size_along_three_directions_for_cubic_lattice(N1, N2, N3) +hamiltonian = guan.hamiltonian_of_finite_size_system_along_three_directions_for_cubic_lattice(N1, N2, N3) LDOS = guan.local_density_of_states_for_cubic_lattice(fermi_energy, hamiltonian, N1=N1, N2=N2, N3=N3) print('cubic lattice:\n', LDOS, '\n') -h00 = guan.finite_size_along_two_directions_for_square_lattice(N2, N3) +h00 = guan.hamiltonian_of_finite_size_system_along_two_directions_for_square_lattice(N2, N3) h01 = np.identity(N2*N3) LDOS = guan.local_density_of_states_for_cubic_lattice_using_dyson_equation(fermi_energy, h00, h01, N3=N3, N2=N2, N1=N1) print(LDOS) \ No newline at end of file diff --git a/Tutorial/calculate_conductance_and_scattering_matrix.py b/Tutorial/06_calculate_conductance_and_scattering_matrix.py similarity index 65% rename from Tutorial/calculate_conductance_and_scattering_matrix.py rename to Tutorial/06_calculate_conductance_and_scattering_matrix.py index 8b10452..0ee5a77 100644 --- a/Tutorial/calculate_conductance_and_scattering_matrix.py +++ b/Tutorial/06_calculate_conductance_and_scattering_matrix.py @@ -2,14 +2,12 @@ import guan import numpy as np # calculate conductance -fermi_energy_array = np.linspace(-5, 5, 400) -h00 = guan.finite_size_along_one_direction(4) +fermi_energy_array = np.linspace(-4, 4, 400) +h00 = guan.hamiltonian_of_finite_size_system_along_one_direction(4) h01 = np.identity(4) conductance_array = guan.calculate_conductance_with_fermi_energy_array(fermi_energy_array, h00, h01) -guan.plot(fermi_energy_array, conductance_array, xlabel='E', ylabel='Conductance', type='-o') +guan.plot(fermi_energy_array, conductance_array, xlabel='E', ylabel='Conductance', type='-') # calculate scattering matrix fermi_energy = 0 -h00 = guan.finite_size_along_one_direction(4) -h01 = np.identity(4) guan.print_or_write_scattering_matrix(fermi_energy, h00, h01) \ No newline at end of file diff --git a/Tutorial/07_calculate_Chern_number_and_Wilson_loop.py b/Tutorial/07_calculate_Chern_number_and_Wilson_loop.py new file mode 100644 index 0000000..d785e0e --- /dev/null +++ b/Tutorial/07_calculate_Chern_number_and_Wilson_loop.py @@ -0,0 +1,13 @@ +import guan +import numpy as np +from math import * + +# calculate Chern number +chern_number = guan.calculate_chern_number_for_square_lattice(guan.hamiltonian_of_one_QAH_model, precision=100) +print('\nChern number=', chern_number) + +# calculate Wilson loop +wilson_loop_array = guan.calculate_wilson_loop(guan.hamiltonian_of_ssh_model) +print('Wilson loop =', wilson_loop_array) +p = np.log(wilson_loop_array)/2/pi/1j +print('p =', p, '\n') \ No newline at end of file diff --git a/Tutorial/read_and_write.py b/Tutorial/08_read_and_write.py similarity index 100% rename from Tutorial/read_and_write.py rename to Tutorial/08_read_and_write.py diff --git a/Tutorial/Fourier_transform_and_calculate_band_structures.py.py b/Tutorial/Fourier_transform_and_calculate_band_structures.py.py deleted file mode 100644 index dc77536..0000000 --- a/Tutorial/Fourier_transform_and_calculate_band_structures.py.py +++ /dev/null @@ -1,9 +0,0 @@ -import guan -import numpy as np -import functools - -# Fourier transform / calculate band structures / plot figures -x_array = np.linspace(-np.pi, np.pi, 100) -hamiltonian_function = functools.partial(guan.one_dimensional_fourier_transform, unit_cell=0, hopping=1) -eigenvalue_array = guan.calculate_eigenvalue_with_one_parameter(x_array, hamiltonian_function) -guan.plot(x_array, eigenvalue_array, xlabel='k', ylabel='E', type='-k') \ No newline at end of file diff --git a/Tutorial/Hamiltonian_of_finite_size_systems.py b/Tutorial/Hamiltonian_of_finite_size_systems.py deleted file mode 100644 index b716715..0000000 --- a/Tutorial/Hamiltonian_of_finite_size_systems.py +++ /dev/null @@ -1,6 +0,0 @@ -import guan - -# Hamiltonian of finite size -print(guan.finite_size_along_one_direction(3), '\n') -print(guan.finite_size_along_two_directions_for_square_lattice(2, 2), '\n') -print(guan.finite_size_along_three_directions_for_cubic_lattice(2, 2, 2), '\n') diff --git a/Tutorial/calculate_Chern_number_and_Wilson_loop.py b/Tutorial/calculate_Chern_number_and_Wilson_loop.py deleted file mode 100644 index 17752a4..0000000 --- a/Tutorial/calculate_Chern_number_and_Wilson_loop.py +++ /dev/null @@ -1,24 +0,0 @@ -import guan -import numpy as np -from math import * - -# calculate Chern number -def hamiltonian_function(kx, ky): # one QAH model with chern number 2 - t1 = 1.0 - t2 = 1.0 - t3 = 0.5 - m = -1.0 - hamiltonian = np.zeros((2, 2), dtype=complex) - hamiltonian[0, 1] = 2*t1*cos(kx)-1j*2*t1*cos(ky) - hamiltonian[1, 0] = 2*t1*cos(kx)+1j*2*t1*cos(ky) - hamiltonian[0, 0] = m+2*t3*sin(kx)+2*t3*sin(ky)+2*t2*cos(kx+ky) - hamiltonian[1, 1] = -(m+2*t3*sin(kx)+2*t3*sin(ky)+2*t2*cos(kx+ky)) - return hamiltonian -chern_number = guan.calculate_chern_number_for_square_lattice(hamiltonian_function, precision=100) -print('Chern number=', chern_number) - -# calculate Wilson loop -wilson_loop_array = guan.calculate_wilson_loop(guan.hamiltonian_of_ssh_model) -print('Wilson loop =', wilson_loop_array) -p = np.log(wilson_loop_array)/2/pi/1j -print('p =', p, '\n') \ No newline at end of file diff --git a/Tutorial/some_models_in_the_reciprocal_space.py b/Tutorial/some_models_in_the_reciprocal_space.py deleted file mode 100644 index 4429b41..0000000 --- a/Tutorial/some_models_in_the_reciprocal_space.py +++ /dev/null @@ -1,9 +0,0 @@ -import guan -import numpy as np - -# Hamiltonian of models in the reciprocal space / calculate band structures / plot figures -x_array = np.linspace(-np.pi, np.pi, 100) -eigenvalue_array = guan.calculate_eigenvalue_with_one_parameter(x_array, guan.hamiltonian_of_square_lattice_in_quasi_one_dimension) -guan.plot(x_array, eigenvalue_array, xlabel='k', ylabel='E', type='-k') -eigenvalue_array = guan.calculate_eigenvalue_with_one_parameter(x_array, guan.hamiltonian_of_graphene_with_zigzag_in_quasi_one_dimension) -guan.plot(x_array, eigenvalue_array, xlabel='k', ylabel='E', type='-k') \ No newline at end of file diff --git a/Tutorial/test_and_Pauli_matrix.py b/Tutorial/test_and_Pauli_matrix.py deleted file mode 100644 index df46229..0000000 --- a/Tutorial/test_and_Pauli_matrix.py +++ /dev/null @@ -1,12 +0,0 @@ -import guan - -# test -print('test') -guan.test() - -# Pauli matrix -print('Pauli matrix') -print('sigma_0:\n', guan.sigma_0(), '\n') -print('sigma_x:\n', guan.sigma_x(), '\n') -print('sigma_y:\n', guan.sigma_y(), '\n') -print('sigma_z:\n', guan.sigma_z(), '\n')