はじめに

分子の世界を探求してみませんか？分子動力学（MD）シミュレーションは、コンピュータを使って分子の動きを予測し、生命現象の解明や新薬開発に役立つ強力なツールです。

この記事では、OpenMMというPythonライブラリを用いて、MDシミュレーションを始めるための入門記事です。
OpenMMは、高速で柔軟性が高く、GPUにも対応しているため、複雑な分子系のシミュレーションも手軽に行えます。

今回は、DNA二重鎖のMDシミュレーションを例に、コードの各部分を丁寧に解説します。
この記事を読めば、あなたも分子シミュレーションの世界への第一歩を踏み出せるでしょう！

今回のコードでできること

今回のコードは、以下のステップでタンパク質の分子動力学シミュレーションを実行し、その結果をグラフで表示します。

1. PDBファイルの読み込み: タンパク質の構造データ（PDBファイル）を読み込みます。
2. 力場の設定: 分子間の相互作用を定義する力場を設定します。
3. 系の構築: タンパク質に水分子とイオンを加えて、シミュレーションの系を構築します。
4. エネルギー最小化: 系の初期構造のエネルギーを最小化し、安定化させます。
5. 初期化: シミュレーションの初期条件（温度など）を設定します。
6. 分子動力学シミュレーション: 設定した条件下で、分子の運動を計算します。
7. データの出力: エネルギーや温度などの物理量を記録します。
8. 結果の可視化: ポテンシャルエネルギーの時間変化をグラフで表示します。

コードの解説

以下にコード全体と、各部分の詳細な解説を記載します。

from openmm.app import *
from openmm import *
from openmm.unit import *
import mdtraj as md
import numpy as np
import matplotlib.pyplot as plt
import sys
import os
import pandas as pd

# Load the PDB file
pdb = PDBFile('/content/drive/My Drive/Colab_Notebooks/1bna.pdb')

# Set up the force field
forcefield = ForceField('amber14-all.xml', 'amber14/tip3p.xml')

# Create the system
modeller = Modeller(pdb.topology, pdb.positions)
modeller.addHydrogens(forcefield)
modeller.addSolvent(forcefield, model='tip3p', boxSize=Vec3(8, 8, 8)*nanometers, ionicStrength=0.1*molar)

# システム作成
system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME, constraints=HBonds, nonbondedCutoff=1.0*nanometer)

# Set up the integrator
temperature = 300*kelvin
pressure = 1*atmosphere
timestep = 2*femtoseconds
integrator = LangevinIntegrator(temperature, 1/picosecond, timestep)
barostat = MonteCarloBarostat(pressure, temperature)
system.addForce(barostat)

# Create the simulation
simulation = Simulation(modeller.topology, system, integrator)
simulation.context.setPositions(modeller.positions)

# Minimize the energy
print("Minimizing energy...")
simulation.minimizeEnergy()

# Equilibrate the system
print("Equilibrating...")
simulation.context.setVelocitiesToTemperature(temperature)
simulation.step(1000)  # Equilibration steps

# Run the simulation
print("Running production simulation...")
simulation.reporters.append(StateDataReporter('/content/drive/My Drive/Colab_Notebooks/output.log', 1000, step=True, potentialEnergy=True,
                                          temperature=True, density=True))
simulation.reporters.append(PDBReporter('/content/drive/My Drive/Colab_Notebooks/trajectory.pdb', 1000))
simulation.step(5000)  # 1 ns = 1,000,000 steps

# データを読み込む
file_path = 'output.log'  # ファイル名
column_names = ["Step", "Potential Energy (kJ/mole)", "Temperature (K)", "Density (g/mL)"]

# CSVとして読み込む
data = pd.read_csv(file_path, skiprows=1, names=column_names)

# "Step" と "Potential Energy (kJ/mole)" を取り出す
steps = data["Step"]
potential_energy = data["Potential Energy (kJ/mole)"]

# グラフを作成
plt.figure(figsize=(10, 6))
plt.plot(steps, potential_energy, marker='o', linestyle='-', color='b', label='Potential Energy')
plt.title('Potential Energy vs. Step', fontsize=16)
plt.xlabel('Step', fontsize=14)
plt.ylabel('Potential Energy (kJ/mole)', fontsize=14)
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend(fontsize=12)
plt.tight_layout()

# グラフを表示
plt.show()

1. ライブラリのインポート

from openmm.app import *
from openmm import *
from openmm.unit import *
import mdtraj as md
import numpy as np
import matplotlib.pyplot as plt
import sys
import os
import pandas as pd

openmm.app: PDBファイルの読み込みや、力場の設定など、OpenMMのアプリケーション層の機能を提供します。
openmm: OpenMMのコア機能を提供します。シミュレーションの実行や、系の設定などを行います。
openmm.unit: 物理量の単位を扱うための機能を提供します。
mdtraj: 分子動力学シミュレーションの軌跡データ（trajectory）を解析するためのライブラリです。今回は直接使用していませんが、MDシミュレーションでは頻繁に使われます。
numpy: 数値計算ライブラリです。
matplotlib.pyplot: グラフ描画ライブラリです。
sys, os: システムやOSに関する機能を提供します。
pandas: データ解析ライブラリです。ログファイルの読み込みに使用します。

2. PDBファイルの読み込み

# Load the PDB file
pdb = PDBFile('/content/drive/My Drive/Colab_Notebooks/1bna_copy.pdb')

PDBFile(): OpenMMの関数で、PDBファイル（Protein Data Bank形式のファイル）を読み込みます。PDBファイルには、タンパク質などの分子の3次元構造情報が格納されています。
'/content/drive/My Drive/Colab_Notebooks/1bna.pdb': PDBファイルのパスを指定します。ご自身の環境に合わせてパスを修正してください。今回は例として「1bna.pdb」というファイルを使用しています。

3. 力場の設定

# Set up the force field
forcefield = ForceField('amber14-all.xml', 'amber14/tip3p.xml')

ForceField(): OpenMMの関数で、力場を設定します。
'amber14-all.xml', 'amber14/tip3p.xml': 使用する力場ファイルを指定します。
amber14-all.xml: タンパク質、核酸、炭水化物、脂質などの生体分子全般に対応したAmber14力場です。(可変)
amber14/tip3p.xml: 水分子モデルとしてTIP3Pを使用します。TIP3Pは、水分子を3つの点でモデル化するシンプルなモデルで、MDシミュレーションで広く使われています。(可変)

4. 系の構築

# Create the system
modeller = Modeller(pdb.topology, pdb.positions)
modeller.addHydrogens(forcefield)
modeller.addSolvent(forcefield, model='tip3p', boxSize=Vec3(8, 8, 8)*nanometers, ionicStrength=0.1*molar)

Modeller(): OpenMMのクラスで、分子系を構築するためのクラスです。
pdb.topology: PDBファイルから読み込んだトポロジー情報（原子の種類、結合情報など）を指定します。
pdb.positions: PDBファイルから読み込んだ原子の座標情報を指定します。
modeller.addHydrogens(forcefield): 系に水素原子を追加します。PDBファイルには水素原子が含まれていない場合があるため、明示的に追加します。力場に基づいて水素原子の最適な位置が自動的に計算されます。
modeller.addSolvent(forcefield, model='tip3p', boxSize=Vec3(8, 8, 8)nanometers, ionicStrength=0.1molar): 系に溶媒（水分子）を追加します。
model='tip3p': 水分子モデルとしてTIP3Pを指定します。(可変)
boxSize=Vec3(8, 8, 8)nanometers: 水分子を入れる箱のサイズを8nm x 8nm x 8nmに指定します。Vec3(x, y, z)は、3次元ベクトルを表すOpenMMのクラスです。nanometersは単位を表します。(可変)
ionicStrength=0.1*molar: イオン強度を0.1 mol/Lに設定します。これにより、系にNa+とCl-イオンが自動的に追加され、電荷が中性になります。(可変)

5. システム作成

# システム作成
system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME, constraints=HBonds, nonbondedCutoff=1.0*nanometer)

forcefield.createSystem(): 力場とトポロジー情報に基づいて、OpenMMのSystemオブジェクトを作成します。Systemオブジェクトは、シミュレーション対象の分子系全体を表し、力場、粒子、相互作用などが定義されています。
modeller.topology: 系のトポロジー情報を指定します。
nonbondedMethod=PME: 非結合相互作用（静電相互作用、ファンデルワールス相互作用）の計算方法としてPME（Particle Mesh Ewald）法を指定します。PME法は、周期境界条件下での静電相互作用を効率的に計算する手法です。
constraints=HBonds: 結合長制約としてH-結合（水素原子と重原子間の結合）に制約を適用します。これにより、シミュレーションのタイムステップを大きくすることができます。
nonbondedCutoff=1.0nanometer: 非結合相互作用のカットオフ距離を1.0nmに設定します。カットオフ距離よりも遠い相互作用は無視されます。(可変)

6. Integrator（積分器）の設定

# Set up the integrator
temperature = 300*kelvin
pressure = 1*atmosphere
timestep = 2*femtoseconds
integrator = LangevinIntegrator(temperature, 1/picosecond, timestep)
barostat = MonteCarloBarostat(pressure, temperature)
system.addForce(barostat)

temperature = 300*kelvin: シミュレーション温度を300ケルビン（約27℃）に設定します。(可変)
pressure = 1*atmosphere: シミュレーション圧力を1気圧に設定します。(可変)
timestep = 2*femtoseconds: シミュレーションのタイムステップを2フェムト秒（1フェムト秒 = 10^-15秒）に設定します。タイムステップは、シミュレーションの精度と計算時間のバランスを考慮して決定する必要があります。(可変)
integrator = LangevinIntegrator(temperature, 1/picosecond, timestep): ランジュバン動力学法による積分器 (LangevinIntegrator) を作成します。ランジュバン動力学法は、一定温度の条件下でのシミュレーションを行うための手法です。
temperature: シミュレーション温度を指定します。
1/picosecond: 摩擦係数を指定します。
timestep: タイムステップを指定します。
barostat = MonteCarloBarostat(pressure, temperature): モンテカルロバロスタット (MonteCarloBarostat) を作成します。バロスタットは、一定圧力の条件下でのシミュレーションを行うための手法です。モンテカルロバロスタットは、系全体の体積をランダムに変化させることで圧力を制御します。
pressure: シミュレーション圧力を指定します。
temperature: シミュレーション温度を指定します。
system.addForce(barostat): バロスタットをシステムに追加します。

7. Simulationオブジェクトの作成

# Create the simulation
simulation = Simulation(modeller.topology, system, integrator)
simulation.context.setPositions(modeller.positions)

simulation = Simulation(modeller.topology, system, integrator): Simulationオブジェクトを作成します。Simulationオブジェクトは、シミュレーションを実行するための中心となるオブジェクトです。
modeller.topology: 系のトポロジー情報を指定します。
system: Systemオブジェクトを指定します。
integrator: 積分器オブジェクト (integrator) を指定します。
simulation.context.setPositions(modeller.positions): シミュレーションの初期構造として、Modellerで構築した構造の原子座標を設定します。Contextオブジェクトは、シミュレーションの現在の状態（座標、速度、力など）を保持するオブジェクトです。

8. エネルギー最小化

# Minimize the energy
print("Minimizing energy...")
simulation.minimizeEnergy()

simulation.minimizeEnergy(): 系のポテンシャルエネルギーを最小化します。エネルギー最小化は、初期構造に歪みがある場合や、局所的なエネルギー極小点にトラップされている場合に、系を安定化させるために行います。

9. 初期化 (Equilibration)

# Equilibrate the system
print("Equilibrating...")
simulation.context.setVelocitiesToTemperature(temperature)
simulation.step(1000)  # Equilibration steps

print("Equilibrating..."): コンソールに「Equilibrating...」と出力し、初期化（平衡化）ステップを開始することを知らせます。
simulation.context.setVelocitiesToTemperature(temperature): 系の原子速度を、設定した温度 (temperature) に対応する速度分布に従ってランダムに設定します。これにより、系を目的の温度に素早く平衡化させることができます。
simulation.step(1000): シミュレーションを1000ステップ実行します。これは、系を平衡状態に近づけるための初期化ステップ（Equilibration）です。初期化ステップでは、系の温度や圧力を目的の値に安定させることが目的です。(可変)

10. 本番シミュレーション (Production Simulation)

# Run the simulation
print("Running production simulation...")
simulation.reporters.append(StateDataReporter('/content/drive/My Drive/Colab_Notebooks/output.log', 1000, step=True, potentialEnergy=True,
                                             temperature=True, density=True))
simulation.reporters.append(PDBReporter('/content/drive/My Drive/Colab_Notebooks/trajectory.pdb', 1000))
simulation.step(5000)  # 10 ps = 5000 steps

print("Running production simulation..."): コンソールに「Running production simulation...」と出力し、本番シミュレーションを開始することを知らせます。
simulation.reporters.append(StateDataReporter(...)): StateDataReporter をシミュレーションに追加します。StateDataReporter は、シミュレーション中に様々な物理量（エネルギー、温度、密度など）を定期的にログファイル (output.log) に出力するための機能です。
'/content/drive/My Drive/Colab_Notebooks/output.log': ログファイルのパスを指定します。ご自身の環境に合わせてパスを修正してください。
1000: 1000ステップごとにログを出力するように設定します。(可変)
step=True, potentialEnergy=True, temperature=True, density=True: ログに出力する物理量を指定します。ここでは、ステップ数、ポテンシャルエネルギー、温度、密度を出力するように設定しています。
simulation.reporters.append(PDBReporter(...)): PDBReporter をシミュレーションに追加します。PDBReporter は、シミュレーション中の原子座標を定期的にPDBファイル (trajectory.pdb) に出力するための機能です。
'/content/drive/My Drive/Colab_Notebooks/trajectory.pdb': トラジェクトリファイルのパスを指定します。ご自身の環境に合わせてパスを修正してください。
1000: 1000ステップごとに座標を出力するように設定します。(可変)
simulation.step(5000): 本番シミュレーションを5000ステップ実行します。今回のタイムステップは2フェムト秒なので、5000ステップは10ピコ秒（0.01ナノ秒）に相当します。(可変)

11. データ読み込みとグラフ作成

# データを読み込む
file_path = 'output.log'  # ファイル名
column_names = ["Step", "Potential Energy (kJ/mole)", "Temperature (K)", "Density (g/mL)"]

# CSVとして読み込む
data = pd.read_csv(file_path, skiprows=1, names=column_names)

# "Step" と "Potential Energy (kJ/mole)" を取り出す
steps = data["Step"]
potential_energy = data["Potential Energy (kJ/mole)"]

# グラフを作成
plt.figure(figsize=(10, 6))
plt.plot(steps, potential_energy, marker='o', linestyle='-', color='b', label='Potential Energy')
plt.title('Potential Energy vs. Step', fontsize=16)
plt.xlabel('Step', fontsize=14)
plt.ylabel('Potential Energy (kJ/mole)', fontsize=14)
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend(fontsize=12)
plt.tight_layout()

# グラフを表示
plt.show()

file_path = 'output.log': ログファイルのパスを指定します。StateDataReporter で指定したファイル名と同じである必要があります。
column_names = ["Step", "Potential Energy (kJ/mole)", "Temperature (K)", "Density (g/mL)"]: ログファイルの列名を定義します。StateDataReporter で出力した物理量の順序と一致させる必要があります。
data = pd.read_csv(file_path, skiprows=1, names=column_names): pandas の read_csv() 関数を使って、ログファイルをCSV形式で読み込み、DataFrame オブジェクト (data) に格納します。skiprows=1 は、ログファイルの1行目（ヘッダー行）をスキップして読み込むことを意味します。names=column_names は、読み込んだデータに列名を付与します。
steps = data["Step"]: DataFrame オブジェクトから、"Step" 列のデータを steps 変数に抽出します。
potential_energy = data["Potential Energy (kJ/mole)"]: DataFrame オブジェクトから、"Potential Energy (kJ/mole)" 列のデータを potential_energy 変数に抽出します。
plt.figure(figsize=(10, 6)): matplotlib.pyplot の figure() 関数を使って、グラフの描画領域を作成します。figsize=(10, 6) は、グラフのサイズを幅10インチ、高さ6インチに設定します。
plt.plot(steps, potential_energy, marker='o', linestyle='-', color='b', label='Potential Energy'): matplotlib.pyplot の plot() 関数を使って、折れ線グラフを描画します。
steps: x軸のデータとして、ステップ数を指定します。
potential_energy: y軸のデータとして、ポテンシャルエネルギーを指定します。
marker='o': データ点を丸いマーカーで表示するように設定します。
linestyle='-': データ点を実線で結ぶように設定します。
color='b': グラフの色を青色に設定します。
label='Potential Energy': 凡例に表示するラベルを "Potential Energy" に設定します。
plt.title('Potential Energy vs. Step', fontsize=16): グラフのタイトルを設定します。
plt.xlabel('Step', fontsize=14): x軸のラベルを設定します。
plt.ylabel('Potential Energy (kJ/mole)', fontsize=14): y軸のラベルを設定します。
plt.grid(True, linestyle='--', alpha=0.6): グラフにグリッド線を表示します。
plt.legend(fontsize=12): 凡例を表示します。
plt.tight_layout(): グラフのレイアウトを調整し、ラベルなどが重ならないようにします。
plt.show(): グラフを表示します。

コードの実行方法

1. OpenMMと関連ライブラリのインストール:

pip install openmm mdtraj numpy matplotlib pandas

2. PDBファイルの準備:

今回のコードでは、1bna.pdb というPDBファイルを読み込んでいます。
PDB ID: 1BNA のPDBファイルをPDBjなどでダウンロードし](https://www.google.com/search?q=https://www.pdbj.org/pdb/search/structid/1bna)%E3%81%AA%E3%81%A9%E3%81%A7%E3%83%80%E3%82%A6%E3%83%B3%E3%83%AD%E3%83%BC%E3%83%89%E3%81%97)、1bna.pdb というファイル名で保存してください。
コード中のPDBファイルのパス ('/content/drive/My Drive/Colab_Notebooks/1bna.pdb') を、ご自身の環境に合わせて修正してください。

3. Pythonコードの実行:

上記のPythonコードを run_md.py など適当なファイル名で保存し、ターミナルから実行します。

python run_md.py

4. 結果の確認:

シミュレーションが終了すると、output.log (ログファイル)と trajectory.pdb (トラジェクトリファイル)が指定したパスに保存されます。
output.log ファイルを読み込み、ポテンシャルエネルギーのグラフが表示されます。

ちなみに今回のシミュレーションでは、以下のポテンシャルエネルギーの経時変化が示されました。

まとめ

この記事では、OpenMMを使って分子動力学シミュレーションを行う基本的なコードを解説しました。
今回のコードをベースに、様々な分子系や条件でのシミュレーションを試してみてください。

OpenMMは、分子動力学シミュレーションの世界への扉を開く、非常に強力なツールです。
この記事が、あなたのMDシミュレーション探求の一助となれば幸いです。

更なるステップ

シミュレーション時間の延長: simulation.step(5000) のステップ数を増やして、より長い時間のシミュレーションを実行してみましょう。その他にも、記事中に「可変」となっているパラメータはチューニングできるので、是非変えて、実行してみてください。
様々な解析: mdtraj ライブラリを使って、RMSD、RMSF、水素結合数など、様々な物理量を解析してみましょう。
系をカスタマイズ: タンパク質の種類を変えたり、溶媒の種類を変えたり、イオンの種類や濃度を変えたりして、シミュレーション系を色々カスタマイズしてみましょう。
GPUの利用: OpenMMはGPUを利用することで、シミュレーションを大幅に高速化できます。GPU環境を構築して、より大規模なシミュレーションに挑戦してみましょう。

興味のある方は、ぜひ実行してみてください！