分子動力学法

分子動力学法（ぶんしどうりきがくほう、英: molecular dynamics、MD法）は、原子ならびに分子の物理的な動きのコンピューターシミュレーション手法である。原子および分子はある時間の間相互作用することが許され、これによって原子の動的発展の光景が得られる。最も一般的なMD法では、原子および分子のトラクジェクトリは、相互作用する粒子の系についての古典力学におけるニュートンの運動方程式を数値的に解くことによって決定される。この系では粒子間の力およびポテンシャルエネルギーは原子間ポテンシャル（分子力学力場）によって定義される。MD法は元々は1950年代末に理論物理学分野で考え出されたが、今日では主に化学物理学、材料科学、生体分子のモデリングに適用されている。系の静的、動的安定構造や、動的過程（ダイナミクス）を解析する手法。
分子の系は莫大な数の粒子から構成されるため、このような複雑系の性質を解析的に探ることは不可能である。MDシミュレーションは数値的手法を用いることによってこの問題を回避する。しかしながら、長いMDシミュレーションは数学的に悪条件であり、数値積分において累積誤差を生成してしまう。これはアルゴリズムとパラメータの適切な選択によって最小化することができるが、完全に取り除くことはできない。
エルゴード仮説に従う系では、単一の分子動力学シミュレーションの展開は系の巨視的熱力学的性質を決定するために使うことができる。エルゴード系の時間平均はミクロカノニカルアンサンブル（小正準集団）平均に対応する。MDは自然の力をアニメーションすることによって未来を予測する、原子スケールの分子の運動についての理解を可能にする「数による統計力学」や「ニュートン力学のラプラス的視点」とも称されている。
MDシミュレーションでは等温、定圧、等温・定圧、定エネルギー、定積、定ケミカルポテンシャル、グランドカノニカルといった様々なアンサンブル（統計集団）の計算が可能である。また、結合長や位置の固定など様々な拘束条件を付加することもできる。計算対象は、バルク、表面、界面、クラスターなど多様な系を扱える。
MD法で扱える系の規模としては、最大で数億原子からなる系の計算例がある。通常の計算規模は数百から数万原子（分子、粒子）程度である。
通常、ポテンシャル関数は、原子-原子の二体ポテンシャルを組み合わせて表現し、これを計算中に変更しない。そのため化学反応のように、原子間結合の生成・開裂を表現するには、何らかの追加の工夫が必要となる。また、ポテンシャルは経験的・半経験的なパラメータから求められる。
こうしたポテンシャル面の精度の問題を回避するため、ポテンシャル面を電子状態の第一原理計算から求める手法もある。このような方法は、第一原理分子動力学法〔量子（ab initio）分子動力学法〕と呼ばれる。この方法では、ポテンシャル面がより正確なものになるが、扱える原子数は格段に減る（スーパーコンピュータを利用しても、最大で約千個程度）。
また第一原理分子動力学法の多くは、電子状態が常に基底状態であることを前提としているものが多く、電子励起状態や電子状態間の非断熱遷移を含む現象の記述は、こうした手法であってもなお困難である。

歴史

モンテカルロシミュレーションの先行する成功に続いて、1950年代末にアルダーとウェインライトによって、1960年代にラーマンによってそれぞれ独立にMD法が開発された。1957年、アルダーおよびウェインライトは剛体球間の弾性衝突を完全にシミュレーションするためにIBM 704計算機を使用した。1960年、ギブソンらはボルン＝マイヤー型の反発相互作用と凝集面積力を用いることによって固体銅の放射線障害をシミュレーションした。1964年、ラーマンはレナード＝ジョーンズ・ポテンシャルを利用した液体アルゴンの画期的シミュレーションを発表した。自己拡散係数といった系の性質の計算は実験データと遜色がなかった。

応用領域

理論物理学分野で始まったMD法は材料科学において人気を得て、1970年代からは生化学および生物物理学での人気を得ている。MDはX線結晶構造解析あるいはNMR分光法から得られた実験的拘束情報に基づいてタンパク質やその他の高分子の三次元構造を洗練するために頻繁に用いられる。物理学において、MDは薄膜成長やイオン-サブプランテーションといった直接観測することができない原子レベルの現象のダイナミクスを調べるために使われる。また、まだ作成されていないあるいは作成することができないナノテクノロジー装置の物理的性質を調べるためにも使われる。生物物理学および構造生物学では、MD法はリガンドドッキング、脂質二重膜のシミュレーション、ホモロジーモデリング、さらにランダムコイルからポリペプチド鎖の折り畳みをシミュレーションすることによってタンパク質構造をab initioに予測するためにも頻繁に適用されている。

シミュレーション設計の制約

分子動力学シミュレーションの設計は利用可能な計算機能力を考慮しなければならない。計算が合理的な時間で終了できるように、シミュレーションサイズ（n = 粒子の数）、時間ステップ、総シミュレーション時間が選択されなければならない。しかしながら、シミュレーションは調べる自然の過程の時間スケールにとって適切なように十分長くなければならない。シミュレーションから統計的に妥当な結論を得るためには、シミュレーションされる時間は自然の過程の速度論と一致しなければならない。さもなければ、MD法は人間が一歩進むよりも短い時間を観察して人間がどうやって歩くのかについて結論付けるのと同じである。タンパク質およびDNAの動力学に関するほとんどの科学論文はナノ秒（10⁻⁹ s）からマイクロ秒（10⁻⁶ s）のシミュレーションからのデータを用いている。これらのシミュレーションを得るためには、複数CPU日からCPU年が必要である。並列アルゴリズムによって負荷をCPU間で分散することができる。この例としては空間的分解アルゴリズムや力分解アルゴリズムがある。
古典的MDシミュレーションの間、CPUを消費するほとんどのタスクは粒子の内部座標の関数としてのポテンシャルの評価である。このエネルギー評価内で最も計算コストが高いのが非結合（非共有結合）部分である。ランダウのO-記法では、全ての対静電相互作用およびファンデルワールス相互作用があらわに考慮されるとすると、一般的な分子動力学シミュレーションはでスケールする。この計算コストは粒子メッシュエバルト (PME) 法（）、P3M法あるいはより球面カットオフ手法（）といった静電的手法を利用することによって低減することができる。
シミュレーションに必要な総CPU時間に影響を与えるもう一つの要素は、積分時間ステップの大きさである。これはポテンシャルの評価の間の時間の長さである。時間ステップは離散化誤差を避けるのに十分小さいように（すなわち系の最も速い振動周波数よりも小さく）選ばれなければならない。古典的MDの典型的な時間ステップは1フェムト秒（10⁻¹⁵ s）のオーダーである。この値はSHAKE〔最も速い原子（例えば水素）の振動を空間に固定する〕といったアルゴリズムを用いることによって延ばすことができる。複数の時間ステップ法が開発されており、これらによってより遅い長距離力の更新の間隔を延ばすことができる。
溶媒中の分子のシミュレーションでは、露な溶媒（陽溶媒）と露でない溶媒（陰溶媒）のどちらかを選択しなければならない。陽溶媒粒子（TIP3P、SPC/E、SPC-f水モデルなど）は力場によって計算コストを掛けて計算しなければならないのに対して、陰溶媒は平均力手法を用いる。陽溶媒は計算コストが高く、シミュレーション中におよそ10倍を超える粒子を含む必要がある。しかし、陽溶媒の粒度と粘度は溶質分子の特定の性質を再現するために必須である。これは運動力学を再現するために特に重要である。
分子動力学シミュレーションの全ての種類において、シミュレーションの箱の大きさは境界条件アーティファクトを避けるのに十分な程大きくなければならない。境界条件は、端において固定された値を選択することによって、あるいは周期境界条件（シミュレーションの一方がその逆側につながることによてバルク相を模倣する）を採用することによってしばしば扱われる。

正準集団（NVT）

正準集団（カノニカルアンサンブル）では、物質の量（N）、容積（V）、温度（T）が保存される。これは等温分子動力学（constant temperature molecular dynamics、CTMD）と呼ばれることもある。NVTでは、吸熱的過程と発熱的過程のエネルギーはサーモスタット（温度調整器）によって交換される。
MDシミュレーションの境界にエネルギーを加えたり取り除いたりするための様々なサーモスタットアルゴリズムが利用可能であり、カノニカルアンサンブルを近似する。温度を制御するための人気のある手法には、速度リスケーリング（V-rescale）、能勢＝フーバー・サーモスタット、能勢＝フーバー・チェイン、ベレンゼン・サーモスタット、アンダーセン・サーモスタット、ランジュバン動力学がある。ベレンゼン・サーモスタットはフライングアイスキューブ効果を発生する可能性があることに留意すべきである。
これらのアルゴリズムを用いてコンホメーションや速度のカノニカル分布を得るのは簡単ではない。これが系の大きさ、サーモスタットの選択、サーモスタットのパラメータ、時間ステップ、積分器にいかに依存するかは、この分野の多くの論文のテーマとなっている。

等温定圧（NPT）集団

等温定圧集団（アンサンブル）では、物質の量（N）、圧力（P）、温度（T）が保存される。サーモスタットに加えて、バロスタット（圧力調整器）が必要である。NPTアンサンブルは、気温と大気圧に開放されているフラスコを用いた実験室条件に最も密接に対応している。
生物膜のシミュレーションでは、等方性圧力制御は適切ではない。脂質二重膜については、圧力制御は定膜面積（NPAT、Aは面積）あるいは定表面張力γ（NPγT）下で行なわれる。

拡張アンサンブル法

レプリカ交換法は拡張アンサンブルである。これは元々、無秩序なスピン系の遅い動力学を扱うために作られた。並列焼きもどし（parallel tempering）法とも呼ばれる。レプリカ交換MD（REMD）法は、複数の温度で走らせた系の非相互作用レプリカの温度を交換することによって多重極小問題を克服しようと試みている。

MDシミュレーションにおけるポテンシャル

分子動力学シミュレーションはポテンシャル関数（シミュレーション中の粒子の相互作用を決定する項の記述）を必要とする。化学および生物学では通常これは力場と呼ばれ、材料物理学では原子間ポテンシャルと呼ばれる。ポテンシャルは多くの段階の物理学的正確性で定義できる。化学で最も一般的に用いられているものは分子力学法に基づいており、粒子-粒子相互作用の古典的取扱い（構造変化やコンホメーション変化は再現できるが、化学反応は再現できない）を具体化している。
完全な量子力学的記述から古典的ポテンシャルへの簡略化は2つの主要な近似を伴う。1つ目はボルン＝オッペンハイマー近似である。この近似では電子のダイナミクスが非常に速く、核の運動に瞬間的反応すると考えることができる、と述べる。結果として、電子の動きと核の動きは別々に扱うことができる。2つ目の近似は、電子よりもかなり重い核を古典ニュートン動力学に従う点粒子として扱う。古典的分子動力学では、電子の影響は単一のポテンシャルエネルギー表面（通常は基底状態を表す）として近似される。
より細かい詳細が必要な時は、量子力学に基づくポテンシャルが用いられる。また、系の大部分を古典的に扱うが、化学的変換が起こる小さな領域を量子系として扱うハイブリッド古典/量子ポテンシャルも開発されている。

経験的ポテンシャル

化学で用いられる経験的ポテンシャルは力場と呼ばれることが多いのに対して、材料化学分野では原子間ポテンシャルと呼ばれる。
化学におけるほとんどの力場は経験的なものであり、化学結合と関連する結合力、結合角、結合二面角、ファンデルワールス力および静電価と関連する非結合力の和から成る。経験的ポテンシャルはアドホックな機能的近似によって限定的に量子力学的効果を表わす。これらのポテンシャルは原子電荷、原子半径の推定値を反映するファンデルワールスパラメータ、平衡結合長、結合角、結合二面角といった自由なパラメータを含む。これらは、詳細な電子構造（量子化学シミュレーション）あるいは弾性係数、格子パラメータ、分光測定といった経験的な物理的性質に対してフィッティングを行うことで得られる。
非結合性相互作用の非局所的な特性のため、これらは系の全ての粒子間の弱い相互作用を少なくとも含む。その計算は通常、MDシミュレーションの速度のボトルネックである。計算コストを下げるため、力場はシフト打ち切り半径、反応場アルゴリズム、粒子メッシュ・エバルト和、あるいはより新しい粒子-粒子-粒子-メッシュ（P3M）法といった数値的近似を用いる。
化学力場は一般にあらかじめ設定された結合様式を用いる（非経験的動力学法を除く）。したがって、化学力場は化学結合の切断の過程や反応を露にモデル化することができない。一方で、結合次数形式に基づいたもののような物理学におけるポテンシャルの多くは、系の複数の異なる接続や結合の切断を記述することができる。こういったポテンシャルの例としては、炭化水素のためのブレナー・ポテンシャルやそれをC-Si-H系とC-O-H系にさらに発展させたものがある。ReaxFFポテンシャルは、結合次数ポテンシャルと化学力場とを組み合わせた完全な反応力場と見なすことができる。

半経験的ポテンシャル

半経験的ポテンシャルは、量子力学からの行列表示を使用する。しかしながら、行列要素の値は特定の原子軌道の重なりの度合いを見積る経験式によって決定される。次に、この行列は異なる原子軌道の占有率を決定するために対角化され、軌道のエネルギー寄与を決定するために再び経験式が使われる。
強結合ポテンシャルとして知られる半経験的ポテンシャルには様々な種類があり、これらはモデル化される原子によって異なる。

分極可能なポテンシャル

ほとんどの古典的力場は分極率の効果を黙示的に含む（例えば量子化学計算から得られた部分電荷を拡大することによって）。これらの部分電荷は原子の質量に関して固定である。しかし、分子動力学シミュレーションはドルーデ粒子や変動電荷といった異なる手法を用いた誘導双極子の導入によって分極率を明示的にモデル化できる。これによって、局所的な化学的環境に応答する原子間の電荷の動的再分配が可能になる。
長年、分極可能MDシミュレーションは次世代シミュレーションとしてもてはやされてきた。水といった均一な液体については、分極率を含めることによって正確性の向上が達成されてきた。タンパク質についても有望な結果が得られている。しかしながら、シミュレーションにおいて分極率をどのように近似するのが最適化についてはいまだ不確かである。

ab-initio法におけるポテンシャル

古典的分子動力学では、単一のポテンシャルエネルギー表面（通常は基底状態）は力場によって表わされる。これはボルン＝オッペンハイマー近似の結果である。励起状態では、化学反応あるいはより正確な表現が必要な時は、電子の振る舞いを密度汎関数法といった量子力学的手法を用いることによって第一原理から得ることができる。これはab initio分子動力学（AIMD）と呼ばれる。電子の自由度を扱うコストから、このシミュレーションの計算コストは古典的分子動力学よりもかなり高い。これはAIMDがより小さな系あるいはより短い時間に制限されることを意味する。
Ab initio量子力学法は、トラジェクトリ中の配座について必要に応じてその場で系のポテンシャルエネルギーを計算するために使うことができる。この計算は反応座標の近傍で大抵行われる。様々な近似を使うことができるが、これらは経験的当て嵌めではなく理論的考察に基づいている。Ab-initio計算は、電子状態の密度やその他の電子的性質といった、経験的手法からは得ることのできない膨大な情報を与える。Ab-initio法を使用する大きな利点は、共有結合の切断あるいは形成を含む反応を調べる能力である。これらの現象は複数の電子状態に対応する。

ハイブリッドQM/MM法

QM（量子力学的）法は非常に強力である。しかしながら、その計算コストは高い。それに対してMM（古典的あるいは分子力学）法は高速だが、いくつかの制限がある（膨大なパラメータ化の必要性、得られたエネルギー推定値がそれ程正確ではないこと、共有結合が切断/形成する反応のシミュレーションに使うことができないこと、化学的反応に関する正確な詳細を与える能力に限界があること）。QM計算の利点（正確性）とMM計算の利点（速さ）を組み合わせた新たな手法が開発されている。これらの手法は混合あるいはハイブリッド量子力学/分子力学法（ハイブリッドQM/MM法）と呼ばれている。
ハイブリッドQM/MM法の最も重要な利点は速さである。最も分かりやすい場合において古典的分子動力学 (MM) を行うコストはO(n²) と見積られる（nは系中の原子の数）。これは主に静電相互作用項（全ての粒子がその他全ての粒子と相互作用する）のためである。しかしながら、打ち切り半径の使用、周期的対表の更新、粒子-メッシュ・エバルト (PME) 法の派生法によって、このコストをO(n) からO(n²) に減らすることができる。言い換えると、2倍の数の原子の系をシミュレーションすると、2倍から4倍の計算力を要することになる。一方で、最も単純なab-initio計算のコストは典型的にO(n³) あるいはそれ以上を見積られる（制限ハートリー＝フォック計算は～O(n2.7) でスケールすることが示唆されている）。この制限を乗り越えるため、系の小さな部分が量子力学的に取り扱われ（典型的には酵素の活性部位）、残りの系が古典的に取り扱われる。
より洗練された実装では、QM/MM法は量子効果に対して敏感な軽い核（例えば水素）と電子状態の両方を扱うために存在する。これによって、水素の波動関数の生成を行うことができる。この方法論は、水素のトンネリングといった現象を調べるために有用である。QM/MM法が新たな発見をもたらした一つの例は、肝臓のアルコール脱水素酵素におけるヒドリド転移の計算である。この場合、水素原子のトンネリングが重要である（反応速度を決定する）。

操舵分子動力学 (SMD)

操舵分子動力学（Steered molecular dynamics; SMD）シミュレーションでは、望む自由度に沿ってタンパク質の構造を引っ張ることによってその構造を操作するためにタンパク質に力を印加する。これらの実験は原子レベルでのタンパク質における構造変化を明らかにするために用いることができる。SMDは機械的な折り畳み構造のほどけや伸長といった出来事をシミュレーションするためにしばしば用いられている。
SMDには2種類の典型的手順がある。1つは引っ張る速度が一定に保たれるもので、もう1つは印加される力が一定のものである。典型的には、調べる系の部分（例えばタンパク質中の1原子）を調和ポテンシャルによって拘束する。次に特定の原子に一定の速度あるいは一定の力を印加する。シミュレーション中で操作される力、距離、角度を変化させることによって望む反応座標に沿って系を動かすために傘サンプリングが用いられる。傘サンプリングによって、系の配置の全て（高エネルギーと低エネルギーのどちらも）が十分にサンプリングされる。次に、それぞれの配置の自由エネルギー変化を平均力ポテンシャル（PMF）として計算することができる。PMFを計算する人気のある手法は一連の傘サンプリングシミュレーションを解析する重みつきヒストグラム解析法（weighted histogram analysis method; WHAM）である。