分子動力学計算

執筆：藤本 和士（関西大学）

分子動力学計算は、原子一つ一つにかかる力を計算することにより、ニュートンの運動方程式に従って粒子の位置を時間発展させる手法です。つまり、原子一つ一つの軌跡を追跡することができます。得られた軌跡データと統計力学理論を組み合わせることで、自由エネルギー、界面張力、赤外吸収スペクトルなどの実験値と比較できるデータを得ることができます。

 

測定できること

分子運動の軌跡 / 動力学解析 / 自由エネルギー / 相互作用解析 / 破壊シミュレーション



 

原理

分子動力学（MD）計算は、原子一つ一つにかかる力を計算することにより、ニュートンの運動方程式に従って粒子の位置を時間発展させる手法です。今、図１の系のMD計算を行うとしましょう。図1はメルト状態のPMMA10量体の分子が435本入った系のスナップショットです。総原子数は66,120原子です。各原子はお互いに相互作用が働いているので、各原子に力がかかります。それにより原子はニュートンの運動方程式に従って運動します。原子iにかかる力をF_iとすると、ニュートンの運動方程式は、
 

となります。ここで、m_i、a_iは原子iの質量と加速度です。今原子は合計で66,120個ありますので、66,120個の連立微分方程式（(1)式）を解かなければいけません。そのような計算は手計算では不可能ですので、計算機を用いて(1)式の連立微分方程式をMD計算では解いています。この方程式を解くに当たってF_iを知る必要があります。
F_iは次式の通り、考えている系のポテンシャルエネルギー(V)の位置微分から計算することができます。
 

r_iは原子iの座標です。図1の分子に働くポテンシャルエネルギーを考えてみましょう。水素(H)、炭素(C)、酸素(O)の原子がありますが、これらの原子間にはファン・デル・ワールス相互作用が働きます。電気陰性度は原子種毎に異なるので、分子中の原子は分極しており電荷を持ちます。そのため各原子間には電気的相互作用（クーロン相互作用）も働きます。この二つの相互作用を総称して分子間相互作用と呼びます。分子内にも様々な力が働きます。H-CやC-Oなどの化学結合、C-C-Oなどの結合角、C-C-C-Cなどの二面角の相互作用が働きます。これらを総称して分子内相互作用と呼びます。代表的なOPLS-AA力場のポテンシャル関数を次式に示します。
  ここで、右辺の第１項、第２項が分子間相互作用で、それぞれファン・デル・ワールス相互作用とクーロン相互作用です。第３項、第４項、第５項は分子内相互作用で、結合、結合角、二面角に由来します。Nは総原子数、ε_ij、σ_ijは原子i、j間のLneard-Jonesパラメータ、q_iは原子iの電荷を表します。N_b、N_a、N_dは化学結合、結合角、二面角の数を表します。k_b、k_aは結合および結合角のバネ定数、V₁、V₂、V₃、V₄は二面角のパラメータを表しています。これらのパラメータ群を力場と呼びます。

これらのパラメータは量子化学計算から決定したり、様々な実験値に合うように決定します。ただ、研究をするたびに毎回これらのパラメータを決定する必要はありません。長年のMD計算の研究の中で、AMBER、OPLS、CHARMM、COMPASなど多くの力場が提案されてきました。そのため、自身の研究に必要な力場は過去の研究例を参考に選択するのが良いと思います。ただし、これらの力場が絶対的に正いというわけでありません。例えば、化学結合が切れるシミュレーションを行いたい場合、化学結合を調和振動子で表現している(3)式は不適切です。このような場合は、モース型関数を選択すべきでしょう。研究で何を対象に、何を明らかにしたいのかを考えた上で力場を選択する必要があります。時には力場開発から行う必要があることも覚えておいてください。

(3)式のポテンシャルを用いて、(2)式に従って各原子の運動を追跡した場合、系はエネルギー一定のアンサンブル（NVEアンサンブル）を生成します。通常、実験室で行われる実験は温度一定、または温度・圧力一定条件下で行われます。MD計算では、熱浴や圧力浴の自由度も含めて(2)式を解くことにより、温度一定条件や温度・圧力一定条件でのシミュレーションを行うことができます。

図１のMD計算をそのまま実行すると系の周りが真空となります。高分子と真空との界面が生成されるので、表面効果が系に影響します。MD計算を用いた多くの研究はバルクの性質を解明するために行われるため、表面効果は取り除かなければなりません。そのため、一般的なMD計算は周期境界条件（Periodic Boundary Condition: PBC）下で実施されます。PBCとはMDセルを上下左右前奥にレプリカを置くことで真空界面を取り除く計算条件です。論文等で示されるMD計算のスナップショットは図１のようにMDセルだけを通常表示します。真空中に塊があるように見えますが、実際はバルクの計算になっていることに注意しておいてください。
   
 

計算の実行

ここでは具体的なMD計算の実行方法について紹介します。多くの優れたMD計算ソフトが公開されており無料で使えます。代表的なMD計算ソフトを表１に示します。それぞれのソフトウェアは、機能面や計算速度面で特徴があります。高度な計算を行う場合は、研究テーマに合わせて選択してください。
 

ソフトウェア名	入手先
MODYLAS	https://www.modylas.org/
GENESIS	https://www.r-ccs.riken.jp/labs/cbrt/
GROMACS	https://www.gromacs.org/
LAMMPS	https://www.lammps.org/
NAMD	https://www.ks.uiuc.edu/Research/namd/

どのMD計算ソフトを用いたとしても、次の三つのファイルを準備しなければいけません。（ファイルの名前はソフトウェア毎に異なります。）
 

1. 座標ファイル：原子座標、原子速度、MDセルサイズなど
2. トポロジーファイル：力場情報、結合情報など
3. 計算条件ファイル：アンサンブル、時間刻み幅、計算時間、温度、圧力など

これら３つを用意して初めてMD計算が実施可能になります。座標ファイルはその名前の通り、原子の座標が記載されています。それ以外にも原子の速度やMDセルサイズなどの情報も記載されています。ただし、原子の座標だけではどのような分子の計算をしているかわかりません。分子の情報が記載されているのがトポロジーファイルです。１番目の原子はC、２番目の原子はC、３番目の原子はOなどの情報が含まれています。加えて、1番目の原子と２番目の原子が結合しているなどの結合情報も含まれます。MD計算を実行するためには、温度、圧力をそれぞれ300 K、1 barに制御する等の計算条件も必要となります。そのような計算条件を含めたものが計算条件ファイルになります。図２に図１の系の各ファイルの一部を表示しております。このようなファイルを作ることができるソフトウェアも多く公開されています。代表的なものの一つにCHARMM-GUI（https://charmm-gui.org/）があります。まず、こちらのwebサービスと無料のMD計算ソフトを使用して練習されるのをおすすめします。
 

 

さらに勉強するために

表１にあるようにMD計算実行のハードルは低くなっています。しかしながら、原理を理解せずに使用した場合間違った結果を出力しかねません。そのため、少しずつでも構いませんので計算原理を理解することを強くお勧めします。

MD計算をよく理解するためには多くの基礎知識が必要となります。まず、上記で紹介しました通りベースとなる学問は解析（古典）力学です。しかし、多数の原子の運動を扱うためアンサンブルを発生します。そのため統計力学の知識も必要となります。計算対象としてタンパク質、ウィルス、高分子など多岐に渡るため、それぞれの分野の知識も必要となります。さらに解析するためにはコンピュータ（プログラミング）の学習も必要となります。全てをカバーする書籍はありませんがMD計算の専門書等を参考文献として記しておきます。