はじめに
このシリーズは計算化学を題材にした Python プログラミングを紹介した「計算化学者のためのPython入門」の1つです。
特に「計算化学関連の Python プログラム開発に携わる予定」の読者を想定してます。
Python3 は Google Colab を使えばブラウザだけで実行 できます。
本記事の要点
- 計算化学（特に量子化学計算）に関する練習問題を Python でやっつけましょう！
Q1. ファイル処理
Keywords: 引数、ファイル処理、list、関数
計算化学では電子状態計算の結果から分子座標のような情報を抽出して、他の解析に利用することがよくあります。
ここでは、サンプルデータ集 のようなテキストファイルから必要な情報を抜き出すファイル処理の練習をしましょう。
Exercise 1.1: 以下のような convert_xyz2lists.py を作ってみましょう
引数 として xyz ファイルを受け取り、
- 分子座標リスト (list_coord)
- 原子名リスト (list_atom)
- コメント行リスト (list_txt)
のような list に変換する 関数 を convert_xyz2lists.py に実装し、最初の分子座標を xyz 形式で出力しましょう。
つまり、sample.xyz を引数として渡すと
$ python3 convert_xyz2lists.py sample.xyz
First Coordinate:
AtomName0 Atom0-x Atom0-y Atom0-z
...
AtomNameN AtomN-x AtomN-y AtomN-z
のように出力される convert_xyz2lists.py を作成してください。xyz データの例は サンプルデータ集 にあります。
サンプルデータ集 の xyz データを利用した場合、目標とする list は以下のようになります。
list_coord = [[[ 0.000000000000, 0.000000000000, 0.350282000000],
[ 0.000000000000, 0.000000000000, -2.319672000000],
[ 0.000000000000, 0.000000000000, 2.805252000000]],
[[ 0.000000000000, 0.002471483472, -1.125393314970],
[ 0.000000000000, 2.480591394046, 0.547634392158],
[ 0.000000000000, -2.409584630412, 0.520323946128]],
[[ 0.000000000000, -0.000584211972, -1.184759840247],
[-0.000000000000, 1.927896891078, 1.100204193427],
[-0.000000000000, -1.871115768559, 1.057451500330]]]
list_coord = [["Au", "Ag", "Cu"],
["Au", "Ag", "Cu"],
["Au", "Ag", "Cu"]]
list_txt = ["mol-1 linear structure", "mol-2", "mol-3"]
Exercise 1.1 のヒント1
まずはアルゴリズム（解析手順）を考えましょう。例えば、
- コマンドライン（python3 convert_xyz2lists.py sample.xyz）の引数から xyz ファイルのパスを取得
- xyz ファイルを上から 1 行ずつ読む
- 1 つの分子座標を読み終わったらそれぞれの情報を list に格納
- すべて読み終わったら先頭の分子座標を出力
という感じでしょうか。
これらのプログラムを書くためには以下のような知識が必要です。
- コマンドライン引数を取得 (sys.argv)
- ファイルの読み込む (with open など)
- 変数の型 (str, int, float, list など)
- 文字列から list への変換 ("hoge fuga".split())
- list へ情報追加 (list.append("hoge"))
わからなければググって調べてみましょう。
Exercise 1.1 のヒント2
例えばこのようなプログラムになると思います。
import sys
# 関数を定義
def convert_xyz2list(file_path):
list_coord, list_atom, list_txt = [], [], []
# ファイルを開く
with open(file_path) as f: # with 構文を使うと close 忘れを防げる
# 先頭行を読み込む（原子数の行）
line = f.readline()
while line: # line が False に相当するまで読み続ける
# Temporary list を初期化
t_coord, t_atom = [], [] # 分子座標 1 つ分のデータ格納 list
# 原子数を取得
natoms = int(line) # line = "3"
# 次の行を 1 行目を読み込む（コメント行）
line = f.readline() # line = "mol-1 linear structure"
list_txt.append(line) # コメント行は文字列として append すればいい
# 座標を読み込むためのループ
for iatom in range(0, natoms)
line = f.readline() # line = "Au 0.000000000000 0.000000000000 0.350282000000"
# line の "文字列" を [リスト] に分割する
# split_line = ["Au", "0.000000000000", "0.000000000000", "0.350282000000"]
split_line = line.split() # 文字列に対する split 関数
# 1 つの分子の原子名を取得
t_atom.append( split_line リストから原子名の文字列を選択 )
# 1 つの分子の座標を取得
# coord_float = [0.0, 0.0, 0.350282] # float の要素を含む list
coord_float = ( split_line から座標を抽出し, float の list に変換 (複数行になる場合あり) )
t_coord.append(coord_float)
# list に append する
list_coord.append(t_coord)
list_atom.append(t_xyz)
# 次の行を 1 行目を読み込む（原子数の行）
line = f.readline()
# ここまできたら while に戻る
return list_coord, list_atom, list_txt
file_path = ( sys.argv を用いて引数を取得 )
# 自分で定義した関数を実行
l_coord, l_atom, l_txt = convert_xyz2list(file_path)
# 出力して結果を確認
print("First Coordinate:")
# 上の例のように出力を整えるためには工夫が必要（検索キーワード: python3 出力 フォーマット など）
# 個人的にはこちらがおすすめです：https://www.javadrive.jp/python/string/index23.html
natoms = len(l_atom[0])
for iatom in range(0, natoms):
print( (原子名 x座標 y座標 z座標 を出力する) )
Exercise 1.1 のヒント3
ヒント 2 の例では convert_xyz2list 関数を外部ファイルから呼び出すたびに "file_path =..." 以下のプログラムが実行されてしまいます。
以下のようにすることで python ファイルをスクリプトとして単体で実行したときのみ、 if 文以降が実行されるようになります。
import sys
# 関数を定義
def convert_xyz2list(file_path):
(ヒント2と同じ)
return list_coord, list_atom, list_txt
if __name__ == '__main__': # この条件式は python ファイルをスクリプトとして単体で実行した場合、True になる。
file_path = ( sys.argv を用いて引数を取得 )
# 自分で定義した関数を実行
l_coord, l_atom, l_txt = convert_xyz2list(file_path)
「if __name__ == __main__」は結構便利です。忘れてしまっても「python if name main」のような検索キーワードを覚えていれば OK です。
Exercise 1.2: サンプルデータ集 の Hessian を読み込んで Hessian 行列を再構築し、3N × 3N 行列として出力しましょう。
Exercise 1.2 のヒント
闇雲にプログラムを書き始める前にファイル処理の手順（アルゴリズム）を考えましょう。まずは自分が手作業で処理しようと思ったらどうするかを考えましょう。それが（効率的な保証はありませんが）アルゴリズムの 1 つになります。どうすれば Hessian 行列を得ることができるでしょうか？
例えば、
- 全要素が 0 で埋められた 9 × 9 行列を生成し、
- ファイル処理をして下三角部分を埋め、
- 最後に対角項に注意しながら下三角行列を転置して足す
などの操作をすれば Hessian 行列を再構築できます。
Q2. Hessian 行列から振動数を求める
Keywords: 荷重化、対角化、分子振動
Hessian 行列はポテンシャルエネルギー $V$ の核座標 $\boldsymbol{ x }$ に関する二次微分として定義されます。
$$ H_{ij} = \frac{\partial^2 V}{\partial x_i \partial x_j} $$
調和近似におけるポテンシャルエネルギーは
$$ V = \frac{1}{2} \sum_{i,j}^{3N} k_{ij} x_i x_j $$
であり、力の定数 $k_{ii}$ は以下のような関係式で振動数 $\nu_i$ と結ばれます。
$$ \nu_i = \frac{1}{2} \sqrt{\frac{k_{ii}}{\mu}}$$
ここで、$\mu$ は換算質量です。
これまではデカルト座標 $\boldsymbol{ x }$ を用いてきましたが、分子振動を表す基準振動ベクトルを得るためには 荷重 Hessian 行列 に変換する必要があります。原子 $a$ の荷重デカルト座標は
$$ \boldsymbol{ \xi } = \sqrt{m_a} \boldsymbol{ x }$$
で表されます。
プログラムを書く前に以下の 2 つの Exercise に挑戦しましょう。
Exercise 2.1: 上記の Hessian 行列をデカルト座標から荷重座標に変換しましょう。
Exercise 2.1 の答え
Exercise 2.2 のヒント にあります。
Exercise 2.2: Hessian 行列 $\boldsymbol{H}$ → 荷重 Hessian 行列 $\boldsymbol{H}^{\text{MW}}$ に変換するための質量行列 $\boldsymbol{M}$ を求めましょう。例えば、異核二原子分子を例にして考えてみましょう。
$$ \boldsymbol{H}^{\text{MW}} = \boldsymbol{M} \boldsymbol{H} \boldsymbol{M}$$
Exercise 2.2 のヒント
荷重 Hessian は以下のように書き表せます。
$$ \boldsymbol{H}^{\text{MW}}_{ij} = \frac{\partial^2 V}{\partial \xi_i \partial \xi_j} = \frac{\partial^2 V}{\partial \sqrt{m_i} x_i \partial \sqrt{m_j} x_j} $$
ここで $\sqrt{m_i}$ は $x_i$ に対応する質量とします。
以下の Exercise にしたがって python プログラムを書いてみましょう。
Exercise 2.3: Exercise 1.2 で読み込んだ Hessian 行列を荷重 Hessian 行列に変換し、対角化しましょう。その後、固有値から振動数を求めましょう。
Exercise 2.3 のヒント
行列の 対角化 は numpy を使うと便利です。例えば、以下のようにすると正方行列 A が対角化され、固有値 (eigval)、固有ベクトル (eigvec) が得られます。
eigval, eigvec = np.linalg.eig(A)
Exercise 1.2 で Hessian 行列を取得できなかった人へ
以下を自分のプログラムにベタ貼りし、list_hess を Hessian 行列として扱ってください。
（改行などは適宜修正してください）
list_hess = \
[[-1.10011E-03, 0.00000E+00, 0.00000E+00, 5.25932E-04, 0.00000E+00, 0.00000E+00, 5.74180E-04, 0.00000E+00, 0.00000E+00], \
[ 0.00000E+00, -1.10011E-03, 0.00000E+00, 0.00000E+00, 5.25932E-04, 0.00000E+00, 0.00000E+00, 5.74180E-04, 0.00000E+00], \
[ 0.00000E+00, 0.00000E+00, 1.15107E-01, 0.00000E+00, 0.00000E+00, -4.74681E-02, 0.00000E+00, 0.00000E+00, -6.76389E-02], \
[ 5.25932E-04, 0.00000E+00, 0.00000E+00, -2.54690E-04, 0.00000E+00, 0.00000E+00, -2.71242E-04, 0.00000E+00, 0.00000E+00], \
[ 0.00000E+00, 5.25932E-04, 0.00000E+00, 0.00000E+00, -2.54690E-04, 0.00000E+00, 0.00000E+00, -2.71242E-04, 0.00000E+00], \
[ 0.00000E+00, 0.00000E+00, -4.74681E-02, 0.00000E+00, 0.00000E+00, 5.34125E-02, 0.00000E+00, 0.00000E+00, -5.94445E-03], \
[ 5.74180E-04, 0.00000E+00, 0.00000E+00, -2.71242E-04, 0.00000E+00, 0.00000E+00, -3.02938E-04, 0.00000E+00, 0.00000E+00], \
[ 0.00000E+00, 5.74180E-04, 0.00000E+00, 0.00000E+00, -2.71242E-04, 0.00000E+00, 0.00000E+00, -3.02938E-04, 0.00000E+00], \
[ 0.00000E+00, 0.00000E+00, -6.76389E-02, 0.00000E+00, 0.00000E+00, -5.94445E-03, 0.00000E+00, 0.00000E+00, 7.35833E-02],]
Exercise 2.4: サンプルデータ集 では直線形の 3 原子分子から Hessian 行列を求めました。この Hessian 行列から振動数を求め、分子構造との関連を考察しましょう。
Exercise 2.4 のヒント
以下の点について考えてみてください。
- 直線形分子の全自由度の数は？
- 負の固有値（= 虚の振動数）を持つ基準振動モードは何を表す？
- 虚の振動数はいくつある？
また、固有値を振動数に変換する場合、単位に注意する必要があります。
- サンプルデータ集 の Hessian 行列の単位は hartree/bohr^2 である
- 単位変換は 1) すべて原子単位に変換した後、2) 目的の単位系に変換するといった手順で行うのがおすすめです。
Q3. 原子間距離行列を求める
Keywords: 直線距離、ループ高速化
ある分子構造における原子間距離行列（distance matrix）を計算しましょう。
原子 $i$ と原子 $j$ の原子間距離は以下のように定義されます。
$$ d_{ij} = \sqrt{ \sum_{k=1}^3 \left( x_{ik} - x_{jk} \right)^2} $$
これは 3 次元空間の 2 点間の距離と同じです。
Python では numpy の linalg.norm() を使えば簡単に計算することができます。
Exercise 3.1: 任意の xyz ファイルを読み込み、distance matrix を格納した csv ファイルを作成しましょう。
Exercise 3.1 のヒント
ループを使って 2 つの原子を選択し、np.linalg.norm で原子間距離を求め、list に格納しましょう。
csv ファイルは np.savetxt() や pandas.DataFrame に対する to_csv() で作成することができます。
Exercise 3.2: 100原子以上の分子系に対して Exercise 3.1 のスクリプトを実行し、小分子系と計算時間を比較しましょう。
python では以下のようにすれば、実行時間を測ることができます（もう少し簡単に書けるかもしれませんが…）。
import time
import datetime
t1 = time.perf_counter()
(Exercise 3.1 で実装した distance matrix を求める関数)
t2 = time.perf_counter()
t = datetime.timedelta(weeks=0, days=0, hours=0, minutes=0, seconds=(t2 - t1), milliseconds=0, microseconds=0)
print("Running time: ", t)
# Output
>> Running time: 0:00:19.516101
100原子以上の分子系に心当たりが無い場合は、サンプルデータ集からコピペしてください。
Exercise 3.3: 100 原子系の場合、実装によってはかなり時間がかかったと思います。計算速度が速くなるような実装を考えてみてください。
Exercise 3.3 のヒント
Python は for loop が遅いことが知られています。また、numpy は C 言語で書かれているため比較的速度が速いです。
このことを踏まえ、以下の点を改善してみてください。
- 二重ループを回避する
- numpy をうまく活用して距離を求める
- 自分で書いたプログラムのどの部分に時間を要するかを調べて対策する
Exercise 3.3 の回答案
natoms = len(xyz)
npxyz = np.array(xyz)
dmat = np.zeros((natoms, natoms))
for i, iatom in enumerate(npxyz):
jatoms = npxyz[i + 1:]
iatoms = np.repeat([iatom], len(jatoms), axis=0)
d_irow = np.linalg.norm(iatoms - jatoms, axis=1)
dmat[i][i + 1:] = d_irow
# Make a full distance matrix from an upper triangular matrix
dmat = dmat + dmat.T - np.diag(dmat)
Q4. 乱数を用いて金属クラスターを生成する
乱数を使って座標を生成し、単一金属クラスターを作ってみましょう。
乱数はimport randomで生成できます。
Exercise 4.1 のプログラム案
metal_radii = xxxxx # 金属原子の半径
natoms = 10 # クラスターサイズ（原子数）
def put_atom_random(r=10.0):
"""
乱数で原子を配置する関数
原点 (0.0, 0.0, 0.0) を中心として半径 r の球体内部に座標を生成するのがよさそう
Parameters
----------
r: float, default: 10.0
In angstrom
"""
coordinate = xxxxx # 実装する
return coordinate
list_coord = []
for iatom in range(natoms):
icoord = put_atom_random()
list_coord.append(icoord)
# Print metal cluster
"""
xyz 形式になるように list_coord を出力する
"""
リンク
-
Project Euler
- プログラムで解く数学の問題集。プログラミング力と数学力が試される。
- 「Python 練習問題」などでググってお気に入りのサイトを見つけてください。