本記事内では、計算機シミュレーションで出力された座標ファイルから分析を行う際に必要となる、
あるいは知っておくと便利なPythonの機能やモジュールを紹介する。
Pythonの基本的な文法については説明を省いている。
また、本記事はあくまで身内向けの内容であり、あまり丁寧なつくりではない。
分析を行うサンプルファイル
サンプルとして、本記事内では以下のファイルに対して分析を行うことにする。
練習用ディレクトリにコピーしておいて欲しい。
Methane
15
1MET C 1 2.127 1.952 3.113 0.000 0.000 0.000
1MET H 2 2.235 1.963 3.128 0.000 0.000 0.000
1MET H 3 2.092 1.856 3.139 0.000 0.000 0.000
1MET H 4 2.071 2.025 3.164 0.000 0.000 0.000
1MET H 5 2.113 1.968 3.005 0.000 0.000 0.000
2MET C 6 0.191 3.025 1.349 0.000 0.000 0.000
2MET H 7 0.085 3.030 1.352 0.000 0.000 0.000
2MET H 8 0.237 3.074 1.428 0.000 0.000 0.000
2MET H 9 0.231 2.926 1.336 0.000 0.000 0.000
2MET H 10 0.229 3.086 1.267 0.000 0.000 0.000
3MET C 11 2.964 0.815 1.378 0.000 0.000 0.000
3MET H 12 2.906 0.726 1.411 0.000 0.000 0.000
3MET H 13 3.018 0.789 1.285 0.000 0.000 0.000
3MET H 14 2.891 0.893 1.363 0.000 0.000 0.000
3MET H 15 3.020 0.839 1.464 0.000 0.000 0.000
3.50000 3.50000 3.50000
.groファイルの内容を軽く解説しておく。
1行目は分子名、2行目が分子数である。
3行目以降に、各原子に対する情報が記載されている。
左から順に分子番号、原子種、原子の通し番号、現在の座標（x,y,z）、移動前の座標（x,y,z）。
最終行は、分子が入っている箱の大きさ（x,y,z）となっている。
なお、数字の単位は nm である。
ファイルの入出力
sample.groを読み取り、原子種（2列目）が "C" である原子の座標のみを抽出するプログラムを
書いてみることにする。
(1) input()を使用する方法
以下のコードによってやりたいことは実現できる。
file = input("type gro's Path:")
carbon = []
with open(file) as f:
f.readline()
num = f.readline()
for i in range(int(num)):
atom = f.readline().split()
if atom[1] == "C":
carbon.append(atom[3:6])
with open("carbon.dat", mode="w") as f:
for i in range(len(carbon)):
f.write(carbon[i][0] + " " + carbon[i][1] + " " + carbon[i][2] + "\n")
コード内容の解説
with open(ファイルパス) as f:でファイルを指定して開き、fとする。
with ~ asの構文を使わない方法もあるが、こちらの方法であればclose()を忘れる心配がない。
今回はinput()を使用して、コード実行時にファイルを選択できるようにしている。
f.readline()でファイルを1行だけ文字列として読み込む。
読み込んだ文字列を、.split()を使って空白ごとに区切られたリストに変換している。
ループ回数に原子数（2行目）を使うことで使いまわしが可能な作りになっている。
ファイル書き込みの際もwith ~ asの構文が利用できるが、mode="w"と指定しなければならない。
mode="w"は、指定したファイルがなければ新規作成し、あれば上書きする。
なお、ファイル内で改行させたいときには "\n" を改行位置に付ける必要がある。
コードの実行
上記のファイルを作業ディレクトリにコピーし、ターミナルで実行する。
（先頭に付いている "$" はターミナルで実行するということを示す記号であり、打ち込む必要はない。）
$ python3 pick_c.py
実行すると読み込むファイルを尋ねられるので、キーボードで打ち込む。
type gro's Path:./sample.gro
これで、以下のようなファイルが生成されたはずだ。
2.127 1.952 3.113
0.191 3.025 1.349
2.964 0.815 1.378
この方法の問題点
先ほどのコマンドを実行した際に気が付いたかもしれないが、input()で文字列を入力する際には
Tab補完が効かない。入力が面倒なだけではなく、打ち間違いを誘発してしまう。
また、履歴に保存されない。ファイル名を少し変えて実行する場合でも毎度打ち込む必要がある。
(2) sysモジュールを使用する方法
sysモジュールを用い、標準入力と標準出力を取り扱うことで先に述べた問題を解決できる。
import sys
carbon = []
with sys.stdin as f:
f.readline()
num = f.readline()
for i in range(int(num)):
atom = f.readline().split()
if atom[1] == "C":
carbon.append(atom[3:6])
for i in range(len(carbon)):
print(carbon[i][0] + " " + carbon[i][1] + " " + carbon[i][2])
コード内容の解説
sys.stdinを使用することで、ターミナルの標準入力を受け取って扱うことができる。
また、標準出力先をファイルに変更することで、print()を利用して書き込みを行っている。
コードの実行
上記のファイルを作業ディレクトリにコピーし、ターミナルで実行する。
$ python3 pick_c+.py < sample.gro > carbon+.dat
<の後ろに読み込むファイルを、>の後ろに出力先のファイルを指定する。
これで、以下のようなファイルが生成される。
2.127 1.952 3.113
0.191 3.025 1.349
2.964 0.815 1.378
標準入出力を利用することで、ファイル名指定の際にTab補完が利用できるようになった。
また、履歴にファイル名も含めて保存されるようになった。
これでファイル名を少しずつ変更しながらの連続実行も楽になったはずだ。
２点間の距離の計算
carbon+.datには、炭素原子3つ（C1, C2, C3 とする）の座標が記載されている。
このデータをもとに、C1-C2, C1-C3, C2-C3 間の距離を計算して表示したい。
(1) 各成分に対して個別に計算する方法
ファイルから座標を読み取り、座標の差から素直に計算すればよい。
import sys
import math
carbons = []
with sys.stdin as f:
for i in range(3):
carbon = f.readline().split()
for i in range(3):
carbon[i] = float(carbon[i])
carbons.append(carbon)
for i in range(len(carbons)-1):
for j in range(i+1,len(carbons)):
dist = math.sqrt((carbons[i][0]-carbons[j][0])**2 + (carbons[i][1]-carbons[j][1])**2 + (carbons[i][2]-carbons[j][2])**2)
print(f"{i+1}-{j+1} {dist}")
コード内容の解説
今回もsysモジュールを使ってファイルを読み取り、print()で書き込む方法をとる。
ファイルから読み取った時点ではstr型なので、float()を使って実数に直す必要がある。
2点間の距離は $\sqrt{dx^2 + dy^2 + dz^2}$ で算出できる。
このコード内では、平方根を計算するためにmathモジュールのsqrt()を用いている。
最終行のように f"{変数}" としてやると、文字列の中に変数を含ませることができる。
コードの実行
上記のファイルを作業ディレクトリにコピーし、ターミナルで実行する。
$ python3 distance.py < carbon+.dat > distance.dat
これで、以下のようなファイルが生成される。
1-2 2.8303923756256832
1-3 2.2368645466366535
2-3 3.5460499150463183
この方法でも目的は問題なく達成できるが、x, y, z で3回も引き算をするのは面倒。
次の方法を使えば、もっとスマートに計算できる。
(2) numpyを使用する方法
numpyモジュールには、pythonで数学的な処理をする際に便利な関数が数多くある。
import sys
import numpy
carbons = []
with sys.stdin as f:
for i in range(3):
carbon = f.readline().split()
for i in range(3):
carbon[i] = float(carbon[i])
carbons.append(carbon)
for i in range(len(carbons)-1):
for j in range(i+1,len(carbons)):
dist = numpy.sqrt(sum(numpy.square(numpy.array(carbons[i])-numpy.array(carbons[j]))))
print(f"{i+1}-{j+1} {dist}")
コード内容の解説
numpyのarray()を利用すると、リストをN次元配列に変換することができる。
配列に変換すれば、要素同士を演算することができるようになる。
これを用いて、2つのリストの要素同士を引き算したリストを算出している。
その後、numpy.square()で各要素を2乗する。
最後はsum()を用いて要素の総和を求め、numpy.sqrt()で平方根をとる。
コードの実行
上記のファイルを作業ディレクトリにコピーし、ターミナルで実行する。
$ python3 distance+.py < carbon+.dat > distance+.dat
これで、以下のようなファイルが生成される。
1-2 2.8303923756256832
1-3 2.2368645466366535
2-3 3.5460499150463183
もっと簡単に計算するには？
さて、結構簡単になったように見えるが、まだ簡単にできる。
2点間の距離、つまり2点の間にあるベクトルの大きさのことを「(L2)ノルム」という。
numpyには、ノルムを直接的に計算できる機能がある。
import sys
import numpy
carbons = []
with sys.stdin as f:
for i in range(3):
carbon = f.readline().split()
for i in range(3):
carbon[i] = float(carbon[i])
carbons.append(carbon)
for i in range(len(carbons)-1):
for j in range(i+1,len(carbons)):
dist = numpy.linalg.norm(numpy.array(carbons[i])-numpy.array(carbons[j]))
print(f"{i+1}-{j+1} {dist}")
numpy.linalg.norm()の引数にベクトルの成分を与えるとL2ノルムを計算してくれる。
追記
もっと詳しい人ならもっと短縮できるようです。
import sys
import numpy as np
import itertools as it
carbons = np.loadtxt(sys.stdin)
for i,j in it.combinations(range(len(carbons)), 2):
dist = np.linalg.norm(carbons[i] - carbons[j])
print(f"{i+1}-{j+1} {dist}")
コード内容の解説
numpy.loattxt()を使えば、ファイルの内容を直接的にN次元配列として、float型で読み取れる。
最初から配列として読み取ってあるので、numpy.array()を使わなくても要素同士の引き算ができる。
また、intertoolsモジュールのcombinations()を使えば、組み合わせを算出できる。
intertools.combinations(n,m)で、「nからm個取り出す際のすべての組み合わせ」を算出する。
これを利用して、距離を求める原子の組み合わせを決定している。