matplotlib for solid state physics

固体物理学にあらわれるグラフを、 matplotlib で描く

グラフの形を見たいだけなので、定数や係数は適当です。

初歩から学ぶ固体物理学 矢口裕之 (5.15)

図5.4 さまざまな温度におけるフェルミ分布関数

本で、T=100K から 1000K とあるのは、T=10K から 100K の誤りと思います。

図5.5 さまざまな温度におけるボース分布関数

# coding: utf-8
# 初歩から学ぶ固体物理学 矢口裕之 (5.15)
# 図5.4 さまざまな温度におけるフェルミ分布関数
import numpy as np
import matplotlib.pyplot as plt
E=np.arange(0, 600)
for T in 1, 10, 50, 100:
    f = 1 / (np.exp((E - 200)/T) + 1)
    plt.plot(E, f)
plt.show()
# 図5.5 さまざまな温度におけるボース分布関数
E=np.arange(0.001, 10, 0.1)
plt.ylim(0, 1000)
for T in 1, 100, 500, 1000:
    f = 1 / (np.exp((E - 0)/T) - 1)
    plt.plot(E, f)
plt.show()

初歩から学ぶ固体物理学 矢口裕之 (7.10)

図7.7 2種類の原子からなる1次元の格子振動における波数と角振動数の関係

# coding: utf-8
# 初歩から学ぶ固体物理学 矢口裕之 (7.10)
# 図7.7 2種類の原子からなる1次元の格子振動における波数と角振動数の関係
import numpy as np
import matplotlib.pyplot as plt
# a=K=1
Ma = 2
Mb = 1
k=np.arange(-np.pi, np.pi, 0.1)
delta = np.sqrt((1/Ma + 1/Mb)**2 - (4* np.sin(k/2)**2)/(Ma * Mb))
omega0 = np.sqrt((1/Ma + 1/Mb) + delta)
omega1 = np.sqrt((1/Ma + 1/Mb) - delta)
plt.plot(k, omega0)
plt.plot(k, omega1)
plt.show()

初歩から学ぶ固体物理学 矢口裕之 (10.26)

図 10.5 ほとんど自由な電子モデルから求められた電子のエネルギー

# coding: utf-8
# 初歩から学ぶ固体物理学 矢口裕之 (10.26)
# 図 10.5 ほとんど自由な電子モデルから求められた電子のエネルギー
import numpy as np
import matplotlib.pyplot as plt
# 逆格子ベクトル Gm
G=np.int(10)
# 周期ポテンシャルの k=Gm のフーリエ係数
V=np.int(4)
# 波数ベクトル
#k=np.arange(-2*G, G, 0.1)
k=np.arange(-1.2* G, 0.2* G, 0.1)
# h, me, V0 は無視しています。
delta=np.sqrt((k**2 - (k+G)**2)**2 + 4*V**2)
E0=k**2 + (k+G)**2 + delta
E1=k**2 + (k+G)**2 - delta
plt.ylim(0, 150)
plt.plot(k, E0)
plt.plot(k, E1)
plt.show()

トポロジカル絶縁体・超伝導体 野村健太郎 (3.41)

図 3.3 (a) 2バンド模型 (3.39) のエネルギー分散

以下を実行すると、ぐりぐり動くよ。

# coding: utf-8
# トポロジカル絶縁体・超伝導体 野村健太郎 (3.41)
# 図 3.3 (a) 2バンド模型 (3.39) のエネルギー分散
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
m=np.int(-1)
kx, ky = np.mgrid[-np.pi:np.pi:0.1, -np.pi:np.pi:0.1]
# t=r=a=1
R=np.sqrt(np.sin(kx)** 2 + np.sin(ky)**2 +
          (m + (1-np.cos(kx)) + (1-np.cos(ky)))**2)
# ε(k)
e = 10
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(kx, ky, e + R)
ax.plot_surface(kx, ky, e - R)
plt.show()

以上