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