はじめに
今回は
・プログラムを作って理解する統計力学-固有エネルギーから各熱力学量を算出する①-
・プログラムを作って理解する統計力学-固有エネルギーから各熱力学量を算出する②-
・プログラムを作って理解する統計力学-固有エネルギーから各熱力学量を算出する③-
で解説した理論をシミュレーションするプログラミングの紹介をします。各記事を合わせて読んでいただきいろいろ数字を変えて統計力学の理解を深めましょう。
コードの解説
コードの全体
import numpy as np
import pandas as pd
kB = 1.380649 * 10 ** -23
h = 6.62607015 * 10 ** -34
e = 1.60217663 * 10 ** -19
h_ = h / (2 * np.pi)
# Temperature = 300
T_list = np.arange(10, 500, 10)
def En(n):
return (n + 1 / 2) * h_ * 10 ** 14
def exp_En(E, T):
z = np.exp(-E / (kB * T))
return z
n = 100
n_list = np.arange(0, n, 1)
T_n = pd.DataFrame()
sum_z = []
sum_E = []
for Temperature in T_list:
En_list = En(n_list)
z = exp_En(En(n_list), Temperature)
p = exp_En(En_list, Temperature) / sum(z)
E = En_list * p
sum_z.append(sum(z))
sum_E.append(sum(E))
zip_list = list(zip(n_list, En_list, z, p, E))
pd.options.display.float_format = "{:.10e}".format
df = pd.DataFrame(zip_list, columns=["n", "En", "z", "p", "E"])
T_n = pd.concat([T_n, df["E"]], axis=1)
T_n.index = n_list
T_n.columns = T_list
T_n.to_csv("Z.csv")
#zip_list2 = list(zip(T_list, sum_E))
#df2 = pd.DataFrame(zip_list2, columns=["Temperature", "<E>"])
#df2.to_csv("AverageE.csv", index=None)
zip_list3 = list(zip(T_list, sum_z, sum_E))
df3 = pd.DataFrame(zip_list3, columns=["Temperature", "Z", "<E>"])
df3["Helm"] = -kB * df3["Temperature"] * np.log(df3["Z"])
df3["S"] = (df3["<E>"] - df3["Helm"]) / df3["Temperature"]
df3.to_csv("state_quantity.csv", index=None)定数と温度
最初に定数と温度の設定をしています。
kB = 1.380649 * 10 ** -23
h = 6.62607015 * 10 ** -34
e = 1.60217663 * 10 ** -19
h_ = h / (2 * np.pi)
# Temperature = 300
T_list = np.arange(10, 500, 10)kBはボルツマン定数で、hはプランク定数、eは素電荷です。
T_listは計算させたい温度の範囲を指定しています。今回の場合10 ~ 500 [K]を10 [K]ずつ設定しています。ここは好きな値を設定してみてください。
固有エネルギーと確率分布を計算するための関数の定義
ここでは固有エネルギーと確率分布を計算するための関数の定義をします。
def En(n):
return (n + 1 / 2) * h_ * 10 ** 14
def exp_En(E, T):
z = np.exp(-E / (kB * T))
return zpythonでは関数の定義にdefを使います。この解説はこちらの記事を参考にしてください。En(n)は固有エネルギーの関数を定義しています。固有エネルギーは任意に決定できますが、今回は解説のために定義しています。なので、この関数は任意の固有エネルギーに変えてみて、結果を比較してみましょう。
$$ E_n = \left(n+\displaystyle\frac{1}{2}\right)\hbar\omega $$
ωは固有振動数で、今回は10^14です。
次にexp_Enです。こちらの関数は固有エネルギーと温度から
$$ z = \exp\left(-\frac{E_n}{k_BT}\right) $$
を計算しています。こちらの関数は変えないでいろいろ実験してみましょう。
シミュレーション部分
ここから上で設定した数値を本格的にシミュレーションを行っていきます。
sum_z = []
sum_E = []
for Temperature in T_list:
En_list = En(n_list)
z = exp_En(En(n_list), Temperature)
p = exp_En(En_list, Temperature) / sum(z)
E = En_list * p
sum_z.append(sum(z))
sum_E.append(sum(E))forから始まるループがシミュレーション部分です。ループの処理はこちらの記事を参考にしてください。
今回は上で設定したT_listの温度と関数を使って、確率分布pを計算します。ここで分母の分配関数はsum(z)としています。
sum()はnumpy形式のlistの合計を計算しています。分配関数はすべての状態の和でしたので、このような処理を行っています。
確率分布を求めたら、期待値Eを計算します。ここでは(固有エネルギー)×(確率)という重みのみ計算しています。最後に分配関数同様、合算しますので、sum()を行います。
これらの計算が終わったらデータを保存します。ここでは、空っぽのlistにappend()で保存しています。
他のコード
他部分はデータを保存したり、加工したりしています。自由エネルギーとエントロピーを計算するコードだけ記載していきます。
df3["Helm"] = -kB * df3["Temperature"] * np.log(df3["Z"])
df3["S"] = (df3["<E>"] - df3["Helm"]) / df3["Temperature"]こちらはpandasというライブラリを使って計算しています。こちらは別の記事で解説します。
最後に
今回は統計力学の理論をプログラミングし、解説しました。固有エネルギーや温度を変えてみて内部エネルギー、自由エネルギー、エントロピーなどの挙動を確認してみて理解を深めましょう。また、シミュレーションを作ってみて、これまで学んできたコードを実践してみましょう。
2023/9/15 J.A.
コメント