注意
本記事では NumPy を利用した秤量計算について解説しますが、本記事より後に投稿した以下の記事では SymPy を利用した秤量計算について解説しています。SymPy版の方がNumPy版と比較してより一般性の高い場合に適用できますので、もし利用される場合には本記事で紹介するNumPy版よりも以下の記事で紹介しているSymPy版の方を推奨します。
前提
複数種類の試薬（反応物）から生成物を作りたいときに、各試薬をどれくらいの質量で混ぜればいいのかを計算する必要がでてきます。例えば、
生成物
| 組成 | 質量 |
|---|---|
| $ \mathrm{La} \mathrm{O} _{0.5} \mathrm{F} _{0.5} \mathrm{Bi} _{1-x} \mathrm{Pb} _{x} \mathrm{S} _2 $ | $\mathrm{2 \ g}$ |
試薬
| 組成 | 分子量 |
|---|---|
| $ \mathrm{La_2 S_3} $ | 373.99 |
| $ \mathrm{Bi_2 O_3} $ | 465.958 |
| $ \mathrm{Bi_2 S_3} $ | 514.14 |
| $ \mathrm{Pb F_2} $ | 245.2 |
| $ \mathrm{Bi F_3} $ | 265.975 |
| $ \mathrm{Bi} $ | 208.98 |
としたときに、まず、
１．生成物に合う各試薬の分子数の比を求め、
２．分子数あたりの質量である分子量と生成物の質量を使って各試薬の質量を求める。
という2stepで計算します。2番目は簡単ですが、1番目では連立方程式を解く必要がでてきます。
1番目はどういうことかというと、
x_1 \mathrm{La_2 S_3} + x_2 \mathrm{Bi_2 O_3} + \cdots + x_6 \mathrm{Bi} \rightarrow \mathrm{La} \mathrm{O} _{0.5} \mathrm{F} _{0.5} \mathrm{Bi} _{1-x} \mathrm{Pb} _{x} \mathrm{S} _2
が成立するような$x_1, x_2, \cdots, x_6$を求めるということです。
$\mathrm{La}$に注目すると、左辺、右辺それぞれの$\mathrm{La}$の個数は$2 x_1$、$1$であり、これが等しくならなければならないため、$2 x_1 = 1$が成立します。他の試薬についても行うと、
\begin{align}
2 x_1 &= 1 \\
3 x_2 &= 0.5 \\
2 x_4 + 3 x_5 &= 0.5 \\
2 x_2 + 2 x_3 + x_5 + x_6 &= 1 - x \\
x_4 &= x \\
3 x_1 + 3 x_3 &= 2
\end{align}
という連立方程式ができます。これを頑張って手で解くと、
\begin{align}
x_1 &= \frac{1}{2} \\
x_2 &= \frac{1}{6} \\
x_3 &= \frac{1}{6} \\
x_4 &= x \\
x_5 &= \frac{1}{6} - \frac{2}{3} x \\
x_6 &= \frac{1}{6} - \frac{1}{3} x\\
\end{align}
となりました。この試料系で組成中の$x$を変化させるだけならExcelでどうにかできますが、これを一々色々な試料系で手で解くのはかなり面倒ですし、解く過程やExcelに式を書き起こす際にも間違いが生じ得ます。
連立方程式Ax=Bに落とし込む
上の連立方程式を$A \boldsymbol{x} = \boldsymbol{b}$の形に落とし込めれば、numpyなどを使って$A, \boldsymbol{b}$から$\boldsymbol{x}$を簡単に得る（＝連立方程式を解く）ことができます。
import numpy as np
x = np.linalg.solve(A, b)
まずは上の連立方程式を$A \boldsymbol{x} = \boldsymbol{b}$の形に落とし込むことを考えます。」上の連立方程式は、
\begin{align}
2 x_1 + 0 x_2 + 0 x_3 + 0 x_4 + 0 x_5 + 0 x_6 &= 1 \\
0 x_1 + 3 x_2 + 0 x_3 + 0 x_4 + 0 x_5 + 0 x_6 &= 0.5 \\
0 x_1 + 0 x_2 + 0 x_3 + 2 x_4 + 3 x_5 + 0 x_6 &= 0.5 \\
0 x_1 + 2 x_2 + 2 x_3 + 0 x_4 + 1 x_5 + 1 x_6 &= 1 - x \\
0 x_1 + 0 x_2 + 0 x_3 + 1 x_4 + 0 x_5 + 0 x_6 &= x \\
3 x_1 + 0 x_2 + 3 x_3 + 0 x_4 + 0 x_5 + 0 x_6 &= 2 \\
\end{align}
と書き直すことができます。さらに行列とベクトルで書き直せば、
\begin{align}
\begin{pmatrix}
2 & 0 & 0 & 0 & 0 & 0 \\
0 & 3 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 2 & 3 & 0 \\
0 & 2 & 2 & 0 & 1 & 1 \\
0 & 0 & 0 & 1 & 0 & 0 \\
3 & 0 & 3 & 0 & 0 & 0 \\
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3 \\
x_4 \\
x_5 \\
x_6 \\
\end{pmatrix}
=
\begin{pmatrix}
1 \\
0.5 \\
0.5 \\
1 - x \\
x \\
2 \\
\end{pmatrix}
\end{align}
となり、左の行列を$A$、中央のベクトルを$\boldsymbol{x}$、右のベクトルを$\boldsymbol{b}$とおくと、連立方程式$A \boldsymbol{x} = \boldsymbol{b}$となります。
どのようなルールで数字が並んでいるのか
コードの中で$A, \boldsymbol{b}$を作れば$\boldsymbol{x}$を求められることは分かりました。ではコード内で$A, \boldsymbol{b}$を作ろうということになるのですが、そのためにはどのようなルールで$A, \boldsymbol{b}$内に数字が並んでいるかを知る必要があります。元々の連立方程式まで辿ればわかりますが、以下のように単位ベクトルを定義すると、ある組成に含まれる各元素の個数はベクトル1個だけで表現できます。
\begin{align}
\mathrm{La} =
\begin{pmatrix}
1 \\
0 \\
0 \\
0 \\
0 \\
0 \\
\end{pmatrix}
, \mathrm{O} =
\begin{pmatrix}
0 \\
1 \\
0 \\
0 \\
0 \\
0 \\
\end{pmatrix}
, \mathrm{F} =
\begin{pmatrix}
0 \\
0 \\
1 \\
0 \\
0 \\
0 \\
\end{pmatrix}
, \mathrm{Bi} =
\begin{pmatrix}
0 \\
0 \\
0 \\
1 \\
0 \\
0 \\
\end{pmatrix}
, \mathrm{Pb} =
\begin{pmatrix}
0 \\
0 \\
0 \\
0 \\
1 \\
0 \\
\end{pmatrix}
, \mathrm{S} =
\begin{pmatrix}
0 \\
0 \\
0 \\
0 \\
0 \\
1 \\
\end{pmatrix}
\end{align}
例えば、今回の生成物である$ \mathrm{La} \mathrm{O} _{0.5} \mathrm{F} _{0.5} \mathrm{Bi} _{1-x} \mathrm{Pb} _{x} \mathrm{S} _2 $は、
\begin{align}
\mathrm{La} + 0.5 \mathrm{O} + 0.5 \mathrm{F} + (1-x) \mathrm{Bi} + x \mathrm{Pb} + 2 \mathrm{S} =
\begin{pmatrix}
1 \\
0.5 \\
0.5 \\
1 - x \\
x \\
2 \\
\end{pmatrix}
= \boldsymbol{b}
\end{align}
つまり$\boldsymbol{b}$となります。1番目の試薬である$\mathrm{La_2 S_3}$は、
\begin{align}
2 \mathrm{La} + 3 \mathrm{S} =
\begin{pmatrix}
2 \\
0 \\
0 \\
0 \\
0 \\
3 \\
\end{pmatrix}
\end{align}
となります。これは$A$の1列目と一致しています。同様に、$A$の1列目は1番目の試薬、2列目は2番目の試薬、・・・、6列目は6番目の試薬をそれぞれベクトルにしたものに一致します。つまり、
１．生成物に含まれる元素全てを単位ベクトルとして定義し、
２．生成物のベクトルを求める（$=\boldsymbol{b}$）
３．各試薬のベクトルを横に並べて行列にする（$=A$）
という手順で$A, \boldsymbol{b}$を作ることができます。
ここまでのコード
ここまでのコードは以下のようになります。
import numpy as np
def main():
x = 0.08
# 反応物（試薬） 組成
reactant = [
{'La': 2, 'S': 3}, # La2S3
{'Bi': 2, 'O': 3}, # Bi2O3
{'Bi': 2, 'S': 3}, # Bi2S3
{'Pb': 1, 'F': 2}, # PbF2
{'Bi': 1, 'F': 3}, # BiF3
{'Bi': 1} # Bi
]
# 生成物 組成
product = {'La': 1, 'O': 0.5, 'F': 0.5, 'Bi': 1 - x, 'Pb': x, 'S': 2} # LaO0.5F0.5Bi(1-x)Pb(x)S2
vec_size = len(product)
b, A = np.zeros(vec_size), np.zeros((vec_size, vec_size)) # ベクトルb & 行列A
index = {} # ベクトルの定義のインデックス
for i, item in enumerate(product.items()): # ベクトルbの作成
b[i] = item[1]
index[item[0]] = i
print(b)
for j, r_dict in enumerate(reactant): # 行列Aの作成
for item in r_dict.items():
A[index[item[0]], j] = item[1]
print(A)
x = np.linalg.solve(A, b) # 連立方程式 Ax = b の解xを求める
print(x)
if __name__ == '__main__':
main()
解説
反応物（試薬）の定義
# 反応物（試薬） 組成
reactant = [
{'La': 2, 'S': 3}, # La2S3
{'Bi': 2, 'O': 3}, # Bi2O3
{'Bi': 2, 'S': 3}, # Bi2S3
{'Pb': 1, 'F': 2}, # PbF2
{'Bi': 1, 'F': 3}, # BiF3
{'Bi': 1} # Bi
]
このコードでの試薬の定義は、試薬それぞれの組成を辞書型で記し、それをreactantというリストにまとめるという方法を採っています。
辞書型というのは、
a = {'La': 2, 'S': 3}
print(a['La'])
2
のように、定義した文字列に対応する数値などを返してくれる働きを持ちます。
生成物（試薬）の定義
# 生成物 組成
product = {'La': 1, 'O': 0.5, 'F': 0.5, 'Bi': 1 - x, 'Pb': x, 'S': 2} # LaO0.5F0.5Bi(1-x)Pb(x)S2
生成物の方も組成を辞書型で記して定義します。
事前に作りたい生成物の組成と質量を決めていて、それに必要な試薬の組成と分子量がわかっているものとします。x = 0.08とxを定義しているので、1 - xは0.92という数値になります。そのため、最終的な結果にxは登場せず、0.08として計算された数値の結果が出てくることに注意してください。
Aとbの準備
vec_size = len(product)
b, A = np.zeros(vec_size), np.zeros((vec_size, vec_size)) # ベクトルb & 行列A
len(product)でproductの大きさ（長さ）が返されます。今回の場合は'La'、'O'、'F'、'Bi'、'Pb'、'S'の6つを定義したので、6がvec_sizeに代入されます。2行目ではbをvec_size個の要素を持ち、全要素が0であるnumpy配列として定義し、Aをvec_size×vec_sizeで全要素が0である二次元numpy配列として定義しています。
bへの代入
index = {} # ベクトルの定義のインデックス
for i, item in enumerate(product.items()): # ベクトルbの作成
b[i] = item[1]
index[item[0]] = i
print(b)
1行目ではindexという空の辞書型を定義しています。3行目のfor文の中でこのindexに「元素名」：「何番目に1をもつ単位ベクトルなのか」という対応を追加していきますが、このindexを利用するのはAへの代入のときになります。
3行目のfor文では、productの対応が順にitemに代入され、iに何回目のループなのかが代入されます。例として、最初のループでは、i = 0、item = ('La', 1)と代入され、
b[0] = 1
index['La'] = 0
となります。これは、「bの0番目に生成物内のLaの個数を代入」し、また「Laのベクトルは0番目に1を持つ単位ベクトルである」と定義していることになります。
このループが終わった後、bとindexは、
print(b)
print(index)
[1. 0.5 0.5 0.92 0.08 2. ]
{'La': 0, 'O': 1, 'F': 2, 'Bi': 3, 'Pb': 4, 'S': 5}
になります。
Aへの代入
for j, r_dict in enumerate(reactant): # 行列Aの作成
for item in r_dict.items():
A[index[item[0]], j] = item[1]
print(A)
2段階のforループになっています。1段階目のforではreactant内の要素を順にr_dictに、何番目のループなのかをjに代入しています。2段階目のforではr_dict内の対応を順にitemに代入しています。大まかに言えば、1段階目は対象とする試薬とAの列を切り替えるループであり、2段階目は対応する元素について個数をAに代入するループです。例として、最初のループでは、1段階目でj = 0、r_dict = {'La': 2, 'S': 3}と代入された後、2段階目で、item = ('La', 2)が代入されます。item[0]は'La'、index['La']は0なので、
A[index[item[0]], j] → A[index['La'], 0] → A[0, 0]となります。
A[0, 0] = 2
A[5, 0] = 3
ループが終わると、Aに値が代入されています。
print(A)
[[2. 0. 0. 0. 0. 0.]
[0. 3. 0. 0. 0. 0.]
[0. 0. 0. 2. 3. 0.]
[0. 2. 2. 0. 1. 1.]
[0. 0. 0. 1. 0. 0.]
[3. 0. 3. 0. 0. 0.]]
Ax=bを解く
x = np.linalg.solve(A, b) # 連立方程式 Ax = b の解xを求める
print(x)
[0.5 0.16666667 0.16666667 0.08 0.11333333 0.14 ]
質量を求める
連立方程式を解いた後は、各試薬の分子量を使って計算するだけです。具体的には、
１．解いた$\boldsymbol{x}$のそれぞれの要素にそれぞれの分子量を掛ける。
２．１．の和で１．を割ることで、和を$\mathrm{1 \ g}$に規格化する。
３．２．に生成物の目標の質量を掛ける。
最終的なPythonコード
上のコードに質量の計算以外にも色々改変して、以下のコードになりました。
import numpy as np
def main():
x = 0.08
# 反応物（試薬） 組成 分子量
reactant = [
[{'La': 2, 'S': 3}, 373.99], # La2S3
[{'Bi': 2, 'O': 3}, 465.958], # Bi2O3
[{'Bi': 2, 'S': 3}, 514.14], # Bi2S3
[{'Pb': 1, 'F': 2}, 245.2], # PbF2
[{'Bi': 1, 'F': 3}, 265.975], # BiF3
[{'Bi': 1}, 208.98] # Bi
]
# 生成物 組成
product = {'La': 1, 'O': 0.5, 'F': 0.5, 'Bi': 1 - x, 'Pb': x, 'S': 2} # LaO0.5F0.5Bi(1-x)Pb(x)S2
product_mass = 2 # 生成物の質量(g)
reactant_mass = calculate_reactant_mass(reactant, product, product_mass, decimal_place=4)
"""
print_result 結果を表示するか
print_equation Ax=bを表示するか
decimal_place 小数第何位まで表示するか
"""
# 結果 [0.87103717 0.36174473 0.39915065 0.09137284 0.14041231 0.1362823 ]
return
def calculate_reactant_mass(reactant, product, product_mass, print_result=True, print_equation=False, decimal_place=10): # reactant:反応物 product:反応物
num_reactant, num_product = len(reactant), len(product)
if num_reactant != num_product:
raise ArithmeticError(f"試薬の種類数（＝{num_reactant}） と 生成物の元素の種類数（＝{num_product}） が一致しません。")
reactant_dict, reactant_mass = np.array([r[0] for r in reactant]), np.array([r[1] for r in reactant]) # それぞれ辞書型・分子量の部分のみ取り出す
b, A = np.zeros(num_product), np.zeros((num_product, num_reactant)) # ベクトルb & 行列A
index = {} # ベクトルの定義
for i, item in enumerate(product.items()): # ベクトルbの作成
b[i] = item[1]
index[item[0]] = i
for j, r_dict in enumerate(reactant_dict): # 行列Aの作成
for item in r_dict.items():
if item[0] in index:
A[index[item[0]], j] = item[1]
else:
raise ArithmeticError(f"生成物に含まれない元素（＝{item[0]}）が試薬に含まれています。")
A_det = np.linalg.det(A) # 解があるかどうかの確認用にAの行列式を計算
if abs(A_det) < 0.0001:
raise ArithmeticError(f"その試薬の組み合わせでは生成物を作られません。")
x = np.linalg.solve(A, b) # 連立方程式 Ax = b の解xを求める
x_mass = x * reactant_mass
x_mass *= product_mass / x_mass.sum() # 規格化してproduct_massを掛ける
if print_equation: # 連立方程式の表示
print("base vector = ")
print(index)
print("A = ")
print(A)
print("x = ")
print(x)
print("b = ")
print(b)
if print_result: # 結果の表示
print('')
print('-----------------------------------------------------------------------')
print('')
product_str = ""
for item in product.items():
if item[1] == 1:
product_str += item[0]
else:
product_str += item[0] + str(item[1])
print(f'product = {product_str:<8}, product_mass = {product_mass} g')
print('')
print(f'index reactant mass')
for i, r_dict in enumerate(reactant_dict):
reactant_str = ""
for item in r_dict.items():
if item[1] == 1:
reactant_str += item[0]
else:
reactant_str += item[0] + str(item[1])
print(f'{i:>3} {reactant_str:<8} {x_mass[i]:>6.{decimal_place}f} g')
print('')
print('-----------------------------------------------------------------------')
print('')
return x_mass
if __name__ == '__main__':
main()
このコードの特徴
reactantとproductの部分を自由にいじることで、試薬と生成物がどんな組成でどんな元素の数であってもこのコード1つで対応することができ、事前に手計算せずに即座に秤量計算ができます。