diff --git a/academic_codes/2019.12.01_0_Metropolis_sampling/Metropolis_sampling_example_1.py b/academic_codes/2019.12.01_Metropolis_sampling/Metropolis_sampling_example_1.py old mode 100755 new mode 100644 similarity index 97% rename from academic_codes/2019.12.01_0_Metropolis_sampling/Metropolis_sampling_example_1.py rename to academic_codes/2019.12.01_Metropolis_sampling/Metropolis_sampling_example_1.py index 8c4005d..644cdae --- a/academic_codes/2019.12.01_0_Metropolis_sampling/Metropolis_sampling_example_1.py +++ b/academic_codes/2019.12.01_Metropolis_sampling/Metropolis_sampling_example_1.py @@ -1,33 +1,33 @@ -""" -This code is supported by the website: https://www.guanjihuan.com -The newest version of this code is on the web page: https://www.guanjihuan.com/archives/1247 -""" - -import random -import math -import numpy as np -from scipy.stats import norm -import matplotlib.pyplot as plt - - -def function(x): # 期望样品的分布，target distribution function - y = (norm.pdf(x, loc=2, scale=1)+norm.pdf(x, loc=-5, scale=1.5))/2 # loc代表了均值,scale代表标准差 - return y - - -T = 100000 # 取T个样品 -pi = [0 for i in range(T)] # 任意选定一个马尔科夫链初始状态 -for t in np.arange(1, T): - pi_star = np.random.uniform(-10, 10) # a proposed distribution， 例如是均匀分布的，或者是一个依赖pi[t - 1]的分布 - alpha = min(1, (function(pi_star) / function(pi[t - 1]))) # 接收率 - u = random.uniform(0, 1) - if u < alpha: - pi[t] = pi_star # 以alpha的概率接收转移 - else: - pi[t] = pi[t - 1] # 不接收转移 - -pi = np.array(pi) # 转成numpy格式 -print(pi.shape) # 查看抽样样品的维度 -plt.plot(pi, function(pi), '*') # 画出抽样样品期望的分布 # 或用plt.scatter(pi, function(pi)) -plt.hist(pi, bins=100, density=1, facecolor='red', alpha=0.7) # 画出抽样样品的分布 # bins是分布柱子的个数，density是归一化，后面两个参数是管颜色的 -plt.show() +""" +This code is supported by the website: https://www.guanjihuan.com +The newest version of this code is on the web page: https://www.guanjihuan.com/archives/1247 +""" + +import random +import math +import numpy as np +from scipy.stats import norm +import matplotlib.pyplot as plt + + +def function(x): # 期望样品的分布，target distribution function + y = (norm.pdf(x, loc=2, scale=1)+norm.pdf(x, loc=-5, scale=1.5))/2 # loc代表了均值,scale代表标准差 + return y + + +T = 100000 # 取T个样品 +pi = [0 for i in range(T)] # 任意选定一个马尔科夫链初始状态 +for t in np.arange(1, T): + pi_star = np.random.uniform(-10, 10) # a proposed distribution， 例如是均匀分布的，或者是一个依赖pi[t - 1]的分布 + alpha = min(1, (function(pi_star) / function(pi[t - 1]))) # 接收率 + u = random.uniform(0, 1) + if u < alpha: + pi[t] = pi_star # 以alpha的概率接收转移 + else: + pi[t] = pi[t - 1] # 不接收转移 + +pi = np.array(pi) # 转成numpy格式 +print(pi.shape) # 查看抽样样品的维度 +plt.plot(pi, function(pi), '*') # 画出抽样样品期望的分布 # 或用plt.scatter(pi, function(pi)) +plt.hist(pi, bins=100, density=1, facecolor='red', alpha=0.7) # 画出抽样样品的分布 # bins是分布柱子的个数，density是归一化，后面两个参数是管颜色的 +plt.show() diff --git a/academic_codes/2019.12.01_0_Metropolis_sampling/Metropolis_sampling_example_2.m b/academic_codes/2019.12.01_Metropolis_sampling/Metropolis_sampling_example_2.m old mode 100755 new mode 100644 similarity index 96% rename from academic_codes/2019.12.01_0_Metropolis_sampling/Metropolis_sampling_example_2.m rename to academic_codes/2019.12.01_Metropolis_sampling/Metropolis_sampling_example_2.m index fb8cffc..dad09f1 --- a/academic_codes/2019.12.01_0_Metropolis_sampling/Metropolis_sampling_example_2.m +++ b/academic_codes/2019.12.01_Metropolis_sampling/Metropolis_sampling_example_2.m @@ -1,18 +1,18 @@ -% This code is supported by the website: https://www.guanjihuan.com -% The newest version of this code is on the web page: https://www.guanjihuan.com/archives/1247 - -clc;clear all;clf; -s=100000; % 取的样品数 -f=[1,2,3,3,3,3,6,5,4,3,2,1]; % 期望得到样品的分布函数 -d=zeros(1,s); % 初始状态 -x=1; -for i=1:s - y=unidrnd(12); % 1到12随机整数 - alpha=min(1,f(y)/f(x)); % 接收率 - u=rand; - if u 0 else width - 1 # 用到周期性边界条件 - bottom_i = i + 1 if i < (width - 1) else 0 - left_j = j - 1 if j > 0 else height - 1 - right_j = j + 1 if j < (height - 1) else 0 - environment = [[top_i, j], [bottom_i, j], [i, left_j], [i, right_j]] - energy = 0 - if angle == None: - for num_i in range(4): - energy += -np.cos(S[i, j] - S[environment[num_i][0], environment[num_i][1]]) - else: - for num_i in range(4): - energy += -np.cos(angle - S[environment[num_i][0], environment[num_i][1]]) - return energy - - -def calculateAllEnergy(S): # 计算整个体系的能量 - energy = 0 - for i in range(S.shape[0]): - for j in range(S.shape[1]): - energy += getEnergy(i, j, S) - return energy/2 # 作用两次要减半 - - -def plot(S): # 画图 - X, Y = np.meshgrid(np.arange(0, S.shape[0]), np.arange(0, S.shape[0])) - U = np.cos(S) - V = np.sin(S) - plt.figure() - plt.quiver(X, Y, U, V, units='inches') - plt.show() - - -if __name__ == '__main__': - main() - +""" +This code is supported by the website: https://www.guanjihuan.com +The newest version of this code is on the web page: https://www.guanjihuan.com/archives/1249 +""" + +import random +import matplotlib.pyplot as plt +import numpy as np +import copy +import math +import time + + +def main(): + size = 30 # 体系大小 + T = 2 # 温度 + ising = get_one_sample(sizeOfSample=size, temperature=T) # 得到符合玻尔兹曼分布的ising模型 + plot(ising) # 画图 + energy_s = [] + for i in range(size): + for j in range(size): + energy_s0 = getEnergy(i=i, j=j, S=ising) # 获取格点的能量 + energy_s = np.append(energy_s, [energy_s0], axis=0) + plt.hist(energy_s, bins=50, density=1, facecolor='red', alpha=0.7) # 画出格点能量分布 # bins是分布柱子的个数，density是归一化，后面两个参数是管颜色的 + plt.show() + + +def get_one_sample(sizeOfSample, temperature): + S = 2 * np.pi * np.random.random(size=(sizeOfSample, sizeOfSample)) # 随机初始状态，角度是0到2pi + print('体系大小：', S.shape) + initialEnergy = calculateAllEnergy(S) # 计算随机初始状态的能量 + print('系统的初始能量是:', initialEnergy) + newS = np.array(copy.deepcopy(S)) + for i00 in range(1000): # 循环一定次数，得到平衡的抽样分布 + newS = Metropolis(newS, temperature) # Metropolis方法抽样，得到玻尔兹曼分布的样品体系 + newEnergy = calculateAllEnergy(newS) + if np.mod(i00, 100) == 0: + print('循环次数%i, 当前系统能量是:' % (i00+1), newEnergy) + print('循环次数%i, 当前系统能量是:' % (i00 + 1), newEnergy) + return newS + + +def Metropolis(S, T): # S是输入的初始状态， T是温度 + delta_max = 0.5 * np.pi # 角度最大的变化度数，默认是90度，也可以调整为其他 + k = 1 # 玻尔兹曼常数 + for i in range(S.shape[0]): + for j in range(S.shape[0]): + delta = (2 * np.random.random() - 1) * delta_max # 角度变化在-90度到90度之间 + newAngle = S[i, j] + delta # 新角度 + energyBefore = getEnergy(i=i, j=j, S=S, angle=None) # 获取该格点的能量 + energyLater = getEnergy(i=i, j=j, S=S, angle=newAngle) # 获取格点变成新角度时的能量 + alpha = min(1.0, math.exp(-(energyLater - energyBefore)/(k * T))) # 这个接受率对应的是玻尔兹曼分布 + if random.uniform(0, 1) <= alpha: + S[i, j] = newAngle # 接受新状态 + else: + pass # 保持为上一个状态 + return S + + +def getEnergy(i, j, S, angle=None): # 计算(i,j)位置的能量，为周围四个的相互能之和 + width = S.shape[0] + height = S.shape[1] + top_i = i - 1 if i > 0 else width - 1 # 用到周期性边界条件 + bottom_i = i + 1 if i < (width - 1) else 0 + left_j = j - 1 if j > 0 else height - 1 + right_j = j + 1 if j < (height - 1) else 0 + environment = [[top_i, j], [bottom_i, j], [i, left_j], [i, right_j]] + energy = 0 + if angle == None: + for num_i in range(4): + energy += -np.cos(S[i, j] - S[environment[num_i][0], environment[num_i][1]]) + else: + for num_i in range(4): + energy += -np.cos(angle - S[environment[num_i][0], environment[num_i][1]]) + return energy + + +def calculateAllEnergy(S): # 计算整个体系的能量 + energy = 0 + for i in range(S.shape[0]): + for j in range(S.shape[1]): + energy += getEnergy(i, j, S) + return energy/2 # 作用两次要减半 + + +def plot(S): # 画图 + X, Y = np.meshgrid(np.arange(0, S.shape[0]), np.arange(0, S.shape[0])) + U = np.cos(S) + V = np.sin(S) + plt.figure() + plt.quiver(X, Y, U, V, units='inches') + plt.show() + + +if __name__ == '__main__': + main() + diff --git a/academic_codes/2019.12.01_1_simulation_of_Ising_model_with_Monte_Carlo_method/get_the_transition_temperature_by_calculating_magnetic_moments_in_Ising_model.py b/academic_codes/2019.12.01_simulation_of_Ising_model_with_Monte_Carlo_method/get_the_transition_temperature_by_calculating_magnetic_moments_in_Ising_model.py old mode 100755 new mode 100644 similarity index 97% rename from academic_codes/2019.12.01_1_simulation_of_Ising_model_with_Monte_Carlo_method/get_the_transition_temperature_by_calculating_magnetic_moments_in_Ising_model.py rename to academic_codes/2019.12.01_simulation_of_Ising_model_with_Monte_Carlo_method/get_the_transition_temperature_by_calculating_magnetic_moments_in_Ising_model.py index 0f29506..9b8e4bd --- a/academic_codes/2019.12.01_1_simulation_of_Ising_model_with_Monte_Carlo_method/get_the_transition_temperature_by_calculating_magnetic_moments_in_Ising_model.py +++ b/academic_codes/2019.12.01_simulation_of_Ising_model_with_Monte_Carlo_method/get_the_transition_temperature_by_calculating_magnetic_moments_in_Ising_model.py @@ -1,76 +1,76 @@ -""" -This code is supported by the website: https://www.guanjihuan.com -The newest version of this code is on the web page: https://www.guanjihuan.com/archives/1249 -""" - -import random -import matplotlib.pyplot as plt -import numpy as np -import copy -import math -import time - - -def main(): - size = 30 # 体系大小 - for T in np.linspace(0.02, 5, 100): - ising, magnetism = get_one_sample(sizeOfSample=size, temperature=T) - print('温度=', T, ' 磁矩=', magnetism, ' 总能量=', calculateAllEnergy(ising)) - plt.plot(T, magnetism, 'o') - plt.show() - - -def get_one_sample(sizeOfSample, temperature): - newS = np.zeros((sizeOfSample, sizeOfSample)) # 初始状态 - magnetism = 0 - for i00 in range(100): - newS = Metropolis(newS, temperature) - magnetism = magnetism + abs(sum(sum(np.cos(newS))))/newS.shape[0]**2 - magnetism = magnetism/100 - return newS, magnetism - - -def Metropolis(S, T): # S是输入的初始状态， T是温度 - k = 1 # 玻尔兹曼常数 - for i in range(S.shape[0]): - for j in range(S.shape[0]): - newAngle = np.random.randint(-1, 1)*np.pi - energyBefore = getEnergy(i=i, j=j, S=S, angle=None) # 获取该格点的能量 - energyLater = getEnergy(i=i, j=j, S=S, angle=newAngle) # 获取格点变成新角度时的能量 - alpha = min(1.0, math.exp(-(energyLater - energyBefore)/(k * T))) # 这个接受率对应的是玻尔兹曼分布 - if random.uniform(0, 1) <= alpha: - S[i, j] = newAngle # 接受新状态 - else: - pass # 保持为上一个状态 - return S - - -def getEnergy(i, j, S, angle=None): # 计算(i,j)位置的能量，为周围四个的相互能之和 - width = S.shape[0] - height = S.shape[1] - top_i = i - 1 if i > 0 else width - 1 # 用到周期性边界条件 - bottom_i = i + 1 if i < (width - 1) else 0 - left_j = j - 1 if j > 0 else height - 1 - right_j = j + 1 if j < (height - 1) else 0 - environment = [[top_i, j], [bottom_i, j], [i, left_j], [i, right_j]] - energy = 0 - if angle == None: - for num_i in range(4): - energy += -np.cos(S[i, j] - S[environment[num_i][0], environment[num_i][1]]) - else: - for num_i in range(4): - energy += -np.cos(angle - S[environment[num_i][0], environment[num_i][1]]) - return energy - - -def calculateAllEnergy(S): # 计算整个体系的能量 - energy = 0 - for i in range(S.shape[0]): - for j in range(S.shape[1]): - energy += getEnergy(i, j, S) - return energy/2 # 作用两次要减半 - - -if __name__ == '__main__': - main() - +""" +This code is supported by the website: https://www.guanjihuan.com +The newest version of this code is on the web page: https://www.guanjihuan.com/archives/1249 +""" + +import random +import matplotlib.pyplot as plt +import numpy as np +import copy +import math +import time + + +def main(): + size = 30 # 体系大小 + for T in np.linspace(0.02, 5, 100): + ising, magnetism = get_one_sample(sizeOfSample=size, temperature=T) + print('温度=', T, ' 磁矩=', magnetism, ' 总能量=', calculateAllEnergy(ising)) + plt.plot(T, magnetism, 'o') + plt.show() + + +def get_one_sample(sizeOfSample, temperature): + newS = np.zeros((sizeOfSample, sizeOfSample)) # 初始状态 + magnetism = 0 + for i00 in range(100): + newS = Metropolis(newS, temperature) + magnetism = magnetism + abs(sum(sum(np.cos(newS))))/newS.shape[0]**2 + magnetism = magnetism/100 + return newS, magnetism + + +def Metropolis(S, T): # S是输入的初始状态， T是温度 + k = 1 # 玻尔兹曼常数 + for i in range(S.shape[0]): + for j in range(S.shape[0]): + newAngle = np.random.randint(-1, 1)*np.pi + energyBefore = getEnergy(i=i, j=j, S=S, angle=None) # 获取该格点的能量 + energyLater = getEnergy(i=i, j=j, S=S, angle=newAngle) # 获取格点变成新角度时的能量 + alpha = min(1.0, math.exp(-(energyLater - energyBefore)/(k * T))) # 这个接受率对应的是玻尔兹曼分布 + if random.uniform(0, 1) <= alpha: + S[i, j] = newAngle # 接受新状态 + else: + pass # 保持为上一个状态 + return S + + +def getEnergy(i, j, S, angle=None): # 计算(i,j)位置的能量，为周围四个的相互能之和 + width = S.shape[0] + height = S.shape[1] + top_i = i - 1 if i > 0 else width - 1 # 用到周期性边界条件 + bottom_i = i + 1 if i < (width - 1) else 0 + left_j = j - 1 if j > 0 else height - 1 + right_j = j + 1 if j < (height - 1) else 0 + environment = [[top_i, j], [bottom_i, j], [i, left_j], [i, right_j]] + energy = 0 + if angle == None: + for num_i in range(4): + energy += -np.cos(S[i, j] - S[environment[num_i][0], environment[num_i][1]]) + else: + for num_i in range(4): + energy += -np.cos(angle - S[environment[num_i][0], environment[num_i][1]]) + return energy + + +def calculateAllEnergy(S): # 计算整个体系的能量 + energy = 0 + for i in range(S.shape[0]): + for j in range(S.shape[1]): + energy += getEnergy(i, j, S) + return energy/2 # 作用两次要减半 + + +if __name__ == '__main__': + main() + diff --git a/academic_codes/2020.01.03_0_bands_of_qusi_1D_square_lattice/bands_of_qusi_1D_square_lattice.py b/academic_codes/2020.01.03_bands_of_qusi_1D_square_lattice/bands_of_qusi_1D_square_lattice.py old mode 100755 new mode 100644 similarity index 96% rename from academic_codes/2020.01.03_0_bands_of_qusi_1D_square_lattice/bands_of_qusi_1D_square_lattice.py rename to academic_codes/2020.01.03_bands_of_qusi_1D_square_lattice/bands_of_qusi_1D_square_lattice.py index 2202a58..187d432 --- a/academic_codes/2020.01.03_0_bands_of_qusi_1D_square_lattice/bands_of_qusi_1D_square_lattice.py +++ b/academic_codes/2020.01.03_bands_of_qusi_1D_square_lattice/bands_of_qusi_1D_square_lattice.py @@ -1,48 +1,48 @@ -""" -This code is supported by the website: https://www.guanjihuan.com -The newest version of this code is on the web page: https://www.guanjihuan.com/archives/3895 -""" - -import numpy as np -import matplotlib.pyplot as plt -from math import * -import cmath -import functools - - -def hamiltonian(k, N, t): # 准一维方格子哈密顿量 - # 初始化为零矩阵 - h00 = np.zeros((N, N), dtype=complex) - h01 = np.zeros((N, N), dtype=complex) - for i in range(N-1): # 原胞内的跃迁h00 - h00[i, i+1] = t - h00[i+1, i] = t - for i in range(N): # 原胞间的跃迁h01 - h01[i, i] = t - matrix = h00 + h01*cmath.exp(1j*k) + h01.transpose().conj()*cmath.exp(-1j*k) - return matrix - - -def main(): - H_k = functools.partial(hamiltonian, N=10, t=1) - k = np.linspace(-pi, pi, 300) - plot_bands_one_dimension(k, H_k) - - -def plot_bands_one_dimension(k, hamiltonian): - dim = hamiltonian(0).shape[0] - dim_k = k.shape[0] - eigenvalue_k = np.zeros((dim_k, dim)) - i0 = 0 - for k0 in k: - matrix0 = hamiltonian(k0) - eigenvalue, eigenvector = np.linalg.eig(matrix0) - eigenvalue_k[i0, :] = np.sort(np.real(eigenvalue[:])) - i0 += 1 - for dim0 in range(dim): - plt.plot(k, eigenvalue_k[:, dim0], '-k') - plt.show() - - -if __name__ == '__main__': - main() +""" +This code is supported by the website: https://www.guanjihuan.com +The newest version of this code is on the web page: https://www.guanjihuan.com/archives/3895 +""" + +import numpy as np +import matplotlib.pyplot as plt +from math import * +import cmath +import functools + + +def hamiltonian(k, N, t): # 准一维方格子哈密顿量 + # 初始化为零矩阵 + h00 = np.zeros((N, N), dtype=complex) + h01 = np.zeros((N, N), dtype=complex) + for i in range(N-1): # 原胞内的跃迁h00 + h00[i, i+1] = t + h00[i+1, i] = t + for i in range(N): # 原胞间的跃迁h01 + h01[i, i] = t + matrix = h00 + h01*cmath.exp(1j*k) + h01.transpose().conj()*cmath.exp(-1j*k) + return matrix + + +def main(): + H_k = functools.partial(hamiltonian, N=10, t=1) + k = np.linspace(-pi, pi, 300) + plot_bands_one_dimension(k, H_k) + + +def plot_bands_one_dimension(k, hamiltonian): + dim = hamiltonian(0).shape[0] + dim_k = k.shape[0] + eigenvalue_k = np.zeros((dim_k, dim)) + i0 = 0 + for k0 in k: + matrix0 = hamiltonian(k0) + eigenvalue, eigenvector = np.linalg.eig(matrix0) + eigenvalue_k[i0, :] = np.sort(np.real(eigenvalue[:])) + i0 += 1 + for dim0 in range(dim): + plt.plot(k, eigenvalue_k[:, dim0], '-k') + plt.show() + + +if __name__ == '__main__': + main() diff --git a/academic_codes/2020.01.03_1_calculation_of_local_currents/calculation_of_local_currents.py b/academic_codes/2020.01.03_calculation_of_local_currents/calculation_of_local_currents.py old mode 100755 new mode 100644 similarity index 97% rename from academic_codes/2020.01.03_1_calculation_of_local_currents/calculation_of_local_currents.py rename to academic_codes/2020.01.03_calculation_of_local_currents/calculation_of_local_currents.py index 75d0fb7..19de32c --- a/academic_codes/2020.01.03_1_calculation_of_local_currents/calculation_of_local_currents.py +++ b/academic_codes/2020.01.03_calculation_of_local_currents/calculation_of_local_currents.py @@ -1,129 +1,129 @@ -""" -This code is supported by the website: https://www.guanjihuan.com -The newest version of this code is on the web page: https://www.guanjihuan.com/archives/3888 -""" - -import numpy as np -import matplotlib.pyplot as plt -import time - - -def matrix_00(width=10): # 电极元胞内跃迁，width不赋值时默认为10 - h00 = np.zeros((width, width)) - for width0 in range(width-1): - h00[width0, width0+1] = 1 - h00[width0+1, width0] = 1 - return h00 - - -def matrix_01(width=10): # 电极元胞间跃迁，width不赋值时默认为10 - h01 = np.identity(width) - return h01 - - -def matrix_LC(width=10, length=300): # 左电极跳到中心区 - h_LC = np.zeros((width, width*length)) - for width0 in range(width): - h_LC[width0, width0] = 1 - return h_LC - - -def matrix_CR(width=10, length=300): # 中心区跳到右电极 - h_CR = np.zeros((width*length, width)) - for width0 in range(width): - h_CR[width*(length-1)+width0, width0] = 1 - return h_CR - - -def matrix_center(width=10, length=300): # 中心区哈密顿量 - hamiltonian = np.zeros((width*length, width*length)) - for length0 in range(length-1): - for width0 in range(width): - hamiltonian[width*length0+width0, width*(length0+1)+width0] = 1 # 长度方向跃迁 - hamiltonian[width*(length0+1)+width0, width*length0+width0] = 1 - for length0 in range(length): - for width0 in range(width-1): - hamiltonian[width*length0+width0, width*length0+width0+1] = 1 # 宽度方向跃迁 - hamiltonian[width*length0+width0+1, width*length0+width0] = 1 - # 中间加势垒 - for j0 in range(6): - for i0 in range(6): - hamiltonian[width*(int(length/2)-3+j0)+int(width/2)-3+i0, width*(int(length/2)-3+j0)+int(width/2)-3+i0]= 1e8 - return hamiltonian - - -def main(): - start_time = time.time() - fermi_energy = 0 - width = 60 - length = 100 - h00 = matrix_00(width) - h01 = matrix_01(width) - G_n = Green_n(fermi_energy+1j*1e-9, h00, h01, width, length) - # 下面是提取数据并画图 - direction_x = np.zeros((width, length)) - direction_y = np.zeros((width, length)) - for length0 in range(length-1): - for width0 in range(width): - direction_x[width0, length0] = G_n[width*length0+width0, width*(length0+1)+width0] - for length0 in range(length): - for width0 in range(width-1): - direction_y[width0, length0] = G_n[width*length0+width0, width*length0+width0+1] - # print(direction_x) - # print(direction_y) - X, Y = np.meshgrid(range(length), range(width)) - plt.quiver(X, Y, direction_x, direction_y) - plt.show() - end_time = time.time() - print('运行时间=', end_time-start_time) - - - -def transfer_matrix(fermi_energy, h00, h01, dim): # 转移矩阵T。dim是传递矩阵h00和h01的维度 - 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) # np.dot()等效于np.matmul() - 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 # a:b代表 a <= x < b - return transfer # 返回转移矩阵 - - -def green_function_lead(fermi_energy, h00, h01, dim): # 电极的表面格林函数 - transfer = transfer_matrix(fermi_energy, h00, h01, dim) - eigenvalue, eigenvector = np.linalg.eig(transfer) - ind = np.argsort(np.abs(eigenvalue)) - temp = np.zeros((2*dim, 2*dim))*(1+0j) - 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_lead(fermi_energy, h00, h01, width, length): # 电极的自能 - h_LC = matrix_LC(width, length) - h_CR = matrix_CR(width, length) - right_lead_surface, left_lead_surface = green_function_lead(fermi_energy, h00, h01, width) - 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) - return right_self_energy, left_self_energy # 返回右电极的自能和左电极的自能 - - -def Green_n(fermi_energy, h00, h01, width, length): # 计算G_n - right_self_energy, left_self_energy = self_energy_lead(fermi_energy, h00, h01, width, length) - hamiltonian = matrix_center(width, length) - green = np.linalg.inv(fermi_energy*np.identity(width*length)-hamiltonian-left_self_energy-right_self_energy) - right_self_energy = (right_self_energy - right_self_energy.transpose().conj())*1j - left_self_energy = (left_self_energy - left_self_energy.transpose().conj())*1j - G_n = np.imag(np.dot(np.dot(green, left_self_energy), green.transpose().conj())) - return G_n - - -if __name__ == '__main__': +""" +This code is supported by the website: https://www.guanjihuan.com +The newest version of this code is on the web page: https://www.guanjihuan.com/archives/3888 +""" + +import numpy as np +import matplotlib.pyplot as plt +import time + + +def matrix_00(width=10): # 电极元胞内跃迁，width不赋值时默认为10 + h00 = np.zeros((width, width)) + for width0 in range(width-1): + h00[width0, width0+1] = 1 + h00[width0+1, width0] = 1 + return h00 + + +def matrix_01(width=10): # 电极元胞间跃迁，width不赋值时默认为10 + h01 = np.identity(width) + return h01 + + +def matrix_LC(width=10, length=300): # 左电极跳到中心区 + h_LC = np.zeros((width, width*length)) + for width0 in range(width): + h_LC[width0, width0] = 1 + return h_LC + + +def matrix_CR(width=10, length=300): # 中心区跳到右电极 + h_CR = np.zeros((width*length, width)) + for width0 in range(width): + h_CR[width*(length-1)+width0, width0] = 1 + return h_CR + + +def matrix_center(width=10, length=300): # 中心区哈密顿量 + hamiltonian = np.zeros((width*length, width*length)) + for length0 in range(length-1): + for width0 in range(width): + hamiltonian[width*length0+width0, width*(length0+1)+width0] = 1 # 长度方向跃迁 + hamiltonian[width*(length0+1)+width0, width*length0+width0] = 1 + for length0 in range(length): + for width0 in range(width-1): + hamiltonian[width*length0+width0, width*length0+width0+1] = 1 # 宽度方向跃迁 + hamiltonian[width*length0+width0+1, width*length0+width0] = 1 + # 中间加势垒 + for j0 in range(6): + for i0 in range(6): + hamiltonian[width*(int(length/2)-3+j0)+int(width/2)-3+i0, width*(int(length/2)-3+j0)+int(width/2)-3+i0]= 1e8 + return hamiltonian + + +def main(): + start_time = time.time() + fermi_energy = 0 + width = 60 + length = 100 + h00 = matrix_00(width) + h01 = matrix_01(width) + G_n = Green_n(fermi_energy+1j*1e-9, h00, h01, width, length) + # 下面是提取数据并画图 + direction_x = np.zeros((width, length)) + direction_y = np.zeros((width, length)) + for length0 in range(length-1): + for width0 in range(width): + direction_x[width0, length0] = G_n[width*length0+width0, width*(length0+1)+width0] + for length0 in range(length): + for width0 in range(width-1): + direction_y[width0, length0] = G_n[width*length0+width0, width*length0+width0+1] + # print(direction_x) + # print(direction_y) + X, Y = np.meshgrid(range(length), range(width)) + plt.quiver(X, Y, direction_x, direction_y) + plt.show() + end_time = time.time() + print('运行时间=', end_time-start_time) + + + +def transfer_matrix(fermi_energy, h00, h01, dim): # 转移矩阵T。dim是传递矩阵h00和h01的维度 + 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) # np.dot()等效于np.matmul() + 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 # a:b代表 a <= x < b + return transfer # 返回转移矩阵 + + +def green_function_lead(fermi_energy, h00, h01, dim): # 电极的表面格林函数 + transfer = transfer_matrix(fermi_energy, h00, h01, dim) + eigenvalue, eigenvector = np.linalg.eig(transfer) + ind = np.argsort(np.abs(eigenvalue)) + temp = np.zeros((2*dim, 2*dim))*(1+0j) + 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_lead(fermi_energy, h00, h01, width, length): # 电极的自能 + h_LC = matrix_LC(width, length) + h_CR = matrix_CR(width, length) + right_lead_surface, left_lead_surface = green_function_lead(fermi_energy, h00, h01, width) + 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) + return right_self_energy, left_self_energy # 返回右电极的自能和左电极的自能 + + +def Green_n(fermi_energy, h00, h01, width, length): # 计算G_n + right_self_energy, left_self_energy = self_energy_lead(fermi_energy, h00, h01, width, length) + hamiltonian = matrix_center(width, length) + green = np.linalg.inv(fermi_energy*np.identity(width*length)-hamiltonian-left_self_energy-right_self_energy) + right_self_energy = (right_self_energy - right_self_energy.transpose().conj())*1j + left_self_energy = (left_self_energy - left_self_energy.transpose().conj())*1j + G_n = np.imag(np.dot(np.dot(green, left_self_energy), green.transpose().conj())) + return G_n + + +if __name__ == '__main__': main() \ No newline at end of file diff --git a/academic_codes/2020.01.03_1_calculation_of_local_currents/calculation_of_local_currents_with_Kwant.py b/academic_codes/2020.01.03_calculation_of_local_currents/calculation_of_local_currents_with_Kwant.py similarity index 100% rename from academic_codes/2020.01.03_1_calculation_of_local_currents/calculation_of_local_currents_with_Kwant.py rename to academic_codes/2020.01.03_calculation_of_local_currents/calculation_of_local_currents_with_Kwant.py diff --git a/academic_codes/2020.01.03_1_calculation_of_local_currents/calculation_of_local_currents_with_guan.py b/academic_codes/2020.01.03_calculation_of_local_currents/calculation_of_local_currents_with_guan.py similarity index 100% rename from academic_codes/2020.01.03_1_calculation_of_local_currents/calculation_of_local_currents_with_guan.py rename to academic_codes/2020.01.03_calculation_of_local_currents/calculation_of_local_currents_with_guan.py diff --git a/academic_codes/2021.02.08_quantum_transport_in_multi_lead_systems/quantum_transport_in_multi_lead_systems_with_guan.py b/academic_codes/2021.02.08_quantum_transport_in_multi_lead_systems/quantum_transport_in_multi_lead_systems_with_guan.py index 0b135dc..3c1a642 100644 --- a/academic_codes/2021.02.08_quantum_transport_in_multi_lead_systems/quantum_transport_in_multi_lead_systems_with_guan.py +++ b/academic_codes/2021.02.08_quantum_transport_in_multi_lead_systems/quantum_transport_in_multi_lead_systems_with_guan.py @@ -55,13 +55,24 @@ def main(): # lead2 lead3 # lead1(L) lead4(R) # lead6 lead5 - transmission_matrix = guan.calculate_six_terminal_transmission_matrix(fermi_energy, h00_for_lead_4=lead_h00, h01_for_lead_4=lead_h01, h00_for_lead_2=lead_h00, h01_for_lead_2=lead_h01, center_hamiltonian=center_hamiltonian, width=width, length=length, internal_degree=1, moving_step_of_leads=0) - transmission_12_array.append(transmission_matrix[0, 1]) - transmission_13_array.append(transmission_matrix[0, 2]) - transmission_14_array.append(transmission_matrix[0, 3]) - transmission_15_array.append(transmission_matrix[0, 4]) - transmission_16_array.append(transmission_matrix[0, 5]) - transmission_1_all_array.append(transmission_matrix[0, 1]+transmission_matrix[0, 2]+transmission_matrix[0, 3]+transmission_matrix[0, 4]+transmission_matrix[0, 5]) + + transmission_12, transmission_13, transmission_14, transmission_15, transmission_16 = guan.calculate_six_terminal_transmissions_from_lead_1(fermi_energy, h00_for_lead_4=lead_h00, h01_for_lead_4=lead_h01, h00_for_lead_2=lead_h00, h01_for_lead_2=lead_h01, center_hamiltonian=center_hamiltonian, width=width, length=length, internal_degree=1, moving_step_of_leads=0) + transmission_12_array.append(transmission_12) + transmission_13_array.append(transmission_13) + transmission_14_array.append(transmission_14) + transmission_15_array.append(transmission_15) + transmission_16_array.append(transmission_16) + transmission_1_all_array.append(\ + transmission_12+transmission_13+transmission_14+transmission_15+transmission_16) + + # transmission_matrix = guan.calculate_six_terminal_transmission_matrix(fermi_energy, h00_for_lead_4=lead_h00, h01_for_lead_4=lead_h01, h00_for_lead_2=lead_h00, h01_for_lead_2=lead_h01, center_hamiltonian=center_hamiltonian, width=width, length=length, internal_degree=1, moving_step_of_leads=0) + # transmission_12_array.append(transmission_matrix[0, 1]) + # transmission_13_array.append(transmission_matrix[0, 2]) + # transmission_14_array.append(transmission_matrix[0, 3]) + # transmission_15_array.append(transmission_matrix[0, 4]) + # transmission_16_array.append(transmission_matrix[0, 5]) + # transmission_1_all_array.append(transmission_matrix[0, 1]+transmission_matrix[0, 2]+transmission_matrix[0, 3]+transmission_matrix[0, 4]+transmission_matrix[0, 5]) + guan.plot(fermi_energy_array, transmission_12_array, xlabel='Fermi energy', ylabel='Transmission_12') guan.plot(fermi_energy_array, transmission_13_array, xlabel='Fermi energy', ylabel='Transmission_13') guan.plot(fermi_energy_array, transmission_14_array, xlabel='Fermi energy', ylabel='Transmission_14') @@ -69,7 +80,7 @@ def main(): guan.plot(fermi_energy_array, transmission_16_array, xlabel='Fermi energy', ylabel='Transmission_16') guan.plot(fermi_energy_array, transmission_1_all_array, xlabel='Fermi energy', ylabel='Transmission_1_all') end_time = time.time() - print('运行时间=', end_time-start_time) + print('运行时间（分）=', (end_time-start_time)/60) if __name__ == '__main__':