Namd等で計算したMDのトラジェクトリを解析するPythonライブラリ:MDTrajの使い方

MDTraj

PythonのライブラリにMDTrajなるものがあり、トラジェクトリがバイナリで書き出される
GromacsやNamdで原子座標を用いた解析をするのに役立ったためPythonを触ったことがない自分でもかなり役に立ったため、ここに使い方を記す。サポートしている形式も

Wide MD format support, including pdb, xtc, trr, dcd, binpos, netcdf, mdcrd, prmtop, and more.(公式ホームページより)

とかなり幅広い。

MdTrajのインストール

基本的には公式ホームページを見てもらうのが一番わかり易い。

• MDTraj-Installation

自分の場合はanacondaを使い、

$ conda install -c conda-forge mdtraj

とするだけ。
先輩いわくpipよりanacondaのほうが依存性関係の問題を一発で解決してくれるらしい。（伝聞）

使いかた

トラジェクトリの読み込み（Namd:初心者版）

Namdのトラジェクトリである.dcdファイルにはトポロジー(原子名、分子の区切りなど)が含まれていないので、

• トラジェクトリファイル（.dcd）
• トポロジーファイル(.pdb)

の2種を読み込む必要がある。

analy.py

import mdtraj as md
traj = md.load_dcd('output.dcd', top='topology.pdb')
print (traj)

result

<mdtraj.Trajectory with 10000 frames, 6441 atoms at 0x110740a90>

トラジェクトリの読み込み（Namd:中級者版）

公式ホームページのmdtraj.load_dcdのページを見ると

mdtraj.load_dcd(filename, top=None, stride=None, atom_indices=None, frame=None)

と、様々なものが指定できる。一つずつ見ていくと、

• filename ： トラジェクトリのファイル名（必須）

  • もちろんパスを含んで指定してもいい ex : '../../name.pdb'

• top : トポロジーファイル（必須）

  • この場合はVMDのAutopsfによって作成されたautopsf.pdbや初期配置を作るときに使った.pdbなどを使うとよい　　　　　　ex : top='../autopsf.psf'

• stride : 間引いて読み込むときの幅（任意）

  • トラジェクトリが2ps刻みで出力しているなどトラジェクトリファイルがかなり大きいという場合に、これを指定すると間引いて（nフレームごとに）読み込んでくれる。指定しなかったらNoneが入る。当たり前だが整数で指定する。（Namdを使っている人はDCDfragなどを調整してそもそもトラジェクトリのアウトプットを5000フレームごとにしておくなどの工夫をしておこう　　　　　ex:stride=5000

• atom_indices（任意）

  • 系の中でも水分子だけ、など特定の分子座標だけ読み込む場合はこれを指定するとプログラムの高速化になるらしいが、自分は使ったことがない。

• frame（任意）

  • 例えば、123frameのみ読み込みたい!！みたいな需要があるのかわからないが、単一フレームのみ読み込む場合にはこれを指定する。これを指定した場合は上のstrideに値を指定していてもstrideのほうが無視される。

mdtraj.iterload

どうやらトラジェクトリのフレーム数が大きいとき、あるまとまりごとに読み込んで区別できるらしい。
時間経過による変化を追うときには便利…かもしれない

example.py

for chunk in md.iterload('output.dcd',chunk=100 top='topology.pdb')
    print chunk

<mdtraj.Trajectory with 100 frames, 423 atoms at 0x110740a90>
<mdtraj.Trajectory with 100 frames, 423 atoms at 0x110740a90>
<mdtraj.Trajectory with 100 frames, 423 atoms at 0x110740a90>

座標などの指定の仕方

基本的にはmdtraj.Trajectoryのページを参考にすれば大丈夫。ここではごく簡単に説明する。

xyz

xyz : np.ndarray, shape=(n_frames, n_atoms, 3)

• n_frames : トラジェクトリのフレーム数
• n_atoms : トポロジーでの原子番号
• ３ : 軸。それぞれ（x,y,z）に対応。

Pythonでは、配列のインデックスは0から始まることに注意！！
つまり、 (x,y,z)はそれぞれ(0,1,2)に対応している。
n_frameやn_atomsも0始まりに注意！

example.py

print ('x: %s\t y: %s\t z: %s' % (traj.xyz[0,0,0],traj.xyz[0,0,1],traj.xyz[0,0,2]))

result

x: 1.5498047	y: 0.50888556	z: 0.4023079

(0フレーム、0番目の原子のx,y,z座標をそれぞれ出力した)

n_atom

系の原子数を表示

n_residues

系の残基の数（ほぼ分子数と等しい）

n_frames

系のフレーム数

example.py

print ('#How many atoms? :   %s' % traj.n_atoms)
print ('#How many residues? :   %s' % traj.n_residues)
print ('#How many frames? :   %s' % traj.n_frames)

result

# How many atoms? :   23449
# How many residues? :   4563
# How many frames? :   10000

topology.atom

トポロジー中での原子名。原子番号を指定する。

example.py

number = 0
while number < 10:
    print("%s : %s" % ((int(number)),traj.topology.atom(number)))
    number += 1

result

0 : TEOM1-C1
1 : TEOM1-H11
2 : TEOM1-H12
3 : TEOM1-H13
4 : TEOM1-O2
5 : TEOM1-C3
6 : TEOM1-H31
7 : TEOM1-H32
8 : TEOM1-C4
9 : TEOM1-H41

このように、(分子名)-(原子名)で出力される。まずはじめにこのようにして目的の原子の番号を調べると良いだろう。

終わりに

これで座標を用いた計算はできるようになっているが、この他にも動径分布関数やオーダーパラメーターなどMDに欠かせない計算もどうやら簡単にできると思うので、分かり次第ここに追記しようと思う。