【分子動力学】MDのトラジェクトリー解析ライブラリ「MDTraj」はじめの一歩

はじめに

この記事では、pythonで使えるトラジェクトリー解析ライブラリであるMDTrajのインストールや簡単な解析方法を解説します。

前提

Anacondaを導入済み

実行環境

CentOS7

参考

インストール

以下を実行するだけでOK

conda install -c conda-forge mdtraj

トラジェクトリーを読み込む、原子を指定するなど

適当なトラジェクトリーファイルを用意してください。
トラジェクトリーの読み込みは

import mdtraj as md
t = md.load('md.xtc', top='md.gro')

のように行います。

座標データはt.xyzで取得できて、

>>>print(t.xyz.shape)
(101, 1960, 3)

これは全部で101フレーム、1960原子、3次元(x,y,z)の座標のリストになっているということ。

各フレームの時間もこのように取得できます。単位はピコセカンドですね。

#最初の10フレームの時間を取得
>>>print(t.time[0:10])
[ 0. 10. 20. 30. 40. 50. 60. 70. 80. 90.]

また、セルの長さも取得できます。

>>>print(t.unitcell_lengths[-1])
[6.9593244 6.9593244 6.9593244]

トラジェクトリーを解析する

RMSDを計算する

import mdtraj as md
import numpy as np
import matplotlib.pyplot as plt
import csv

# トラジェクトリーの読み込み
t = md.load('md.xtc', top='md.gro')

# 各フレームの時間のリストを取得
time = t.time
# トラジェクトリーの最初のフレームを参照構造としてtのRMSDを計算する
rmsds = md.rmsd(t,t,0)

# 時間とRMSDをCSVファイルに書き込む
with open('rmsd.csv', 'w') as f:
        writer = csv.writer(f)
        for i in range(len(time)):
                writer.writerow([time[i], rmsds[i]])

# 時間vsRMSDでグラフを描画する
plt.plot(t.time, rmsds, label="rmsd")
plt.xlabel('time [ps]')
plt.ylabel('rmsd [nm]')
plt.show()

もしgroファイルを参照構造としたいなら、下記のように書けばよいです。

t = md.load('md.xtc', top='md.gro')
r = md.load('md.gro')

rmsds = md.rmsd(t,r,0)

↓でもう少し突っ込んだ内容も書きました