分子動力学法

分子動力学法（ぶんしどうりきがくほう、英: molecular dynamics、MD法）は、原子ならびに分子の物理的な動きのコンピューターシミュレーション手法である。原子および分子はある時間の間相互作用することが許され、これによって原子の動的発展の光景が得られる。最も一般的なMD法では、原子および分子のトラクジェクトリは、相互作用する粒子の系についての古典力学におけるニュートンの運動方程式を数値的に解くことによって決定される。この系では粒子間の力およびポテンシャルエネルギーは原子間ポテンシャル（分子力学力場）によって定義される。MD法は元々は1950年代末に理論物理学分野で考え出されたが^[1][2]、今日では主に化学物理学、材料科学、生体分子のモデリングに適用されている。系の静的、動的安定構造や、動的過程（ダイナミクス）を解析する手法。

分子の系は莫大な数の粒子から構成されるため、このような複雑系の性質を解析的に探ることは不可能である。MDシミュレーションは 数値的手法を用いることによってこの問題を回避する。しかしながら、長いMDシミュレーションは数学的に悪条件であり、数値積分において累積誤差を生成してしまう。これはアルゴリズムとパラメータの適切な選択によって最小化することができるが、完全に取り除くことはできない。

エルゴード仮説に従う系では、単一の分子動力学シミュレーションの展開は系の巨視的熱力学的性質を決定するために使うことができる。エルゴード系の時間平均はミクロカノニカルアンサンブル（小正準集団）平均に対応する。MDは自然の力をアニメーションすることによって未来を予測する、原子スケールの分子の運動についての理解を可能にする「数による統計力学」や「ニュートン力学のラプラス的視点」とも称されている^[3][4]。

MDシミュレーションでは等温、定圧、等温・定圧、定エネルギー、定積、定ケミカルポテンシャル、グランドカノニカルといった様々なアンサンブル（統計集団）の計算が可能である。また、結合長や位置の固定など様々な拘束条件を付加することもできる。計算対象は、バルク、表面、界面、クラスターなど多様な系を扱える。

MD法で扱える系の規模としては、最大で数億原子からなる系の計算例がある。通常の計算規模は数百から数万原子（分子、粒子）程度である。

通常、ポテンシャル関数は、原子-原子の二体ポテンシャルを組み合わせて表現し、これを計算中に変更しない。そのため化学反応のように、原子間結合の生成・開裂を表現するには、何らかの追加の工夫が必要となる。また、ポテンシャルは経験的・半経験的なパラメータから求められる。

こうしたポテンシャル面の精度の問題を回避するため、ポテンシャル面を電子状態の第一原理計算から求める手法もある。このような方法は、第一原理分子動力学法〔量子（ab initio）分子動力学法〕と呼ばれる。この方法では、ポテンシャル面がより正確なものになるが、扱える原子数は格段に減る（スーパーコンピュータを利用しても、最大で約千個程度）。

また第一原理分子動力学法の多くは、電子状態が常に基底状態であることを前提としているものが多く、電子励起状態や電子状態間の非断熱遷移を含む現象の記述は、こうした手法であってもなお困難である。

歴史[編集]

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

年表[編集]

• 1957年：剛体球の分子動力学法（Alder and Wainwright←最初のMD）
• 1964年：質点系への拡張 (Rahman)
• 1971年：剛体系への拡張 (Rahman and Stillinger)
• 1977年：拘束系への拡張（Rychaert等）
• 1980年：定圧条件の導入（Andersenの方法、Parrinello-Rahman法）
• 1983年：非平衡系への拡張 (Gillan and Dixon)
• 1984年：定温条件の導入（能勢＝フーバーの方法）
• 1985年：第一原理分子動力学法（→カー・パリネロ法）
• 1991年：大正準集団（集合）系への拡張 (Cagin and Pettit)

応用領域[編集]

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

シミュレーション設計の制約[編集]

分子動力学シミュレーションの設計は利用可能な計算機能力を考慮しなければならない。計算が合理的な時間で終了できるように、シミュレーションサイズ（n = 粒子の数）、時間ステップ、総シミュレーション時間が選択されなければならない。しかしながら、シミュレーションは調べる自然の過程の時間スケールにとって適切なように十分長くなければならない。シミュレーションから統計的に妥当な結論を得るためには、シミュレーションされる時間は自然の過程の速度論と一致しなければならない。さもなければ、MD法は人間が一歩進むよりも短い時間を観察して人間がどうやって歩くのかについて結論付けるのと同じである。タンパク質およびDNAの動力学に関するほとんどの科学論文はナノ秒（10⁻⁹ s）からマイクロ秒（10⁻⁶ s）のシミュレーションからのデータを用いている。これらのシミュレーションを得るためには、複数CPU日からCPU年が必要である。並列アルゴリズムによって負荷をCPU間で分散することができる。この例としては空間的分解アルゴリズムや力分解アルゴリズムがある^[6]。

古典的MDシミュレーションの間、CPUを消費するほとんどのタスクは粒子の内部座標の関数としてのポテンシャルの評価である。このエネルギー評価内で最も計算コストが高いのが非結合（非共有結合）部分である。ランダウのO-記法では、全ての対静電相互作用およびファンデルワールス相互作用があらわに考慮されるとすると、一般的な分子動力学シミュレーションは${\displaystyle O(n^{2})}$でスケールする。この計算コストは粒子メッシュエバルト (PME) 法（${\displaystyle O(n\log(n))}$）、P3M法あるいはより球面カットオフ手法（${\displaystyle O(n)}$）といった静電的手法を利用することによって低減することができる。

シミュレーションに必要な総CPU時間に影響を与えるもう一つの要素は、積分時間ステップの大きさである。これはポテンシャルの評価の間の時間の長さである。時間ステップは離散化誤差を避けるのに十分小さいように（すなわち系の最も速い振動周波数よりも小さく）選ばれなければならない。古典的MDの典型的な時間ステップは1フェムト秒（10⁻¹⁵ s）のオーダーである。この値はSHAKE〔最も速い原子（例えば水素）の振動を空間に固定する〕といったアルゴリズムを用いることによって延ばすことができる。複数の時間ステップ法が開発されており、これらによってより遅い長距離力の更新の間隔を延ばすことができる^[7][8][9]。

溶媒中の分子のシミュレーションでは、露な溶媒（陽溶媒）と露でない溶媒（陰溶媒）のどちらかを選択しなければならない。陽溶媒粒子（TIP3P、SPC/E、SPC-f水モデルなど）は力場によって計算コストを掛けて計算しなければならないのに対して、陰溶媒は平均力手法を用いる。陽溶媒は計算コストが高く、シミュレーション中におよそ10倍を超える粒子を含む必要がある。しかし、陽溶媒の粒度と粘度は溶質分子の特定の性質を再現するために必須である。これは運動力学を再現するために特に重要である。

分子動力学シミュレーションの全ての種類において、シミュレーションの箱の大きさは境界条件アーティファクトを避けるのに十分な程大きくなければならない。境界条件は、端において固定された値を選択することによって、あるいは周期境界条件（シミュレーションの一方がその逆側につながることによてバルク相を模倣する）を採用することによってしばしば扱われる。

小正準集団（NVE）[編集]

小正準（ミクロカノニカル、NVE）集団（アンサンブル）において、系はモル（N）、容積（V）、エネルギー（E）の変化から分離される。これは熱交換のない断熱過程に対応する。ミクロカノニカル分子動力学トラクジェクトリは全エネルギーが保存されたポテンシャルエネルギーと運動エネルギーの交換として見ることができる。座標${\displaystyle X}$と${\displaystyle V}$を持つ速度N個の粒子の系では、一次微分方程式の対をニュートンの記法で以下のように書くことができる。

${\displaystyle F(X)=-\nabla U(X)=M{\dot {V}}(t)}$

${\displaystyle V(t)={\dot {X}}(t).}$

系のポテンシャルエネルギー関数${\displaystyle U(X)}$は、粒子の座標${\displaystyle X}$の関数である。これは物理学では「ポテンシャル」、化学では「力場」と単に呼ばれる。最初の方程式はニュートンの法則（系中の個々の粒子に作用する力${\displaystyle F}$は${\displaystyle U(X)}$の負の勾配として計算できる）から来ている。

全ての時間ステップについて、個々の粒子の位置${\displaystyle X}$および速度${\displaystyle V}$はVerlet法といったシンプレティック法を用いて積分することができる。${\displaystyle X}$および${\displaystyle V}$の時間発展はトラジェクトリと呼ばれる。初期位置（例えば理論上の知見から）および初期速度（例えばランダム化ガウシアンから）が与えられれば、未来（あるいは過去）の全ての位置および速度を計算することができる。

よくある混乱の源の一つはMDにおける温度の意味である。一般に、我々が経験しているのは膨大な数の粒子を含む巨視的温度である。しかし温度は統計的量である。もし、十分大きな数の原子が存在すれば、統計的温度は「瞬間温度」から見積ることができる。これは、系の運動エネルギーをnk_BT/2（kBはボルツマン定数、nは系の自由度の数）と同じと見なすことで得られる。

温度に関連した現象はMDシミュレーションで使われる少数の原子が原因で生じる。例えば、500原子を含む基質と100 eVの蒸着エネルギーから開始される銅薄膜の成長のシミュレーションを考える。現実世界では、蒸着した原子からの100 eVは多数の原子（${\displaystyle 10^{10}}$以上）の間ですばやく輸送、共有され、温度に大きな変化は生じない。しかしながら、わずか500原子しかない時は、基質は蒸着によってほぼすぐに蒸発する。生物物理学シミュレーションでも似た事例が起こる。NVEにおける系の温度はタンパク質といった高分子が発熱的なコンホメーション変化や結合を起こす時に自然に上昇する。

正準集団（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）法^[10]は、複数の温度で走らせた系の非相互作用レプリカの温度を交換することによって多重極小問題を克服しようと試みている。

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

分子動力学シミュレーションはポテンシャル関数（シミュレーション中の粒子の相互作用を決定する項の記述）を必要とする。化学および生物学では通常これは力場と呼ばれ、材料物理学では原子間ポテンシャルと呼ばれる。ポテンシャルは多くの段階の物理学的正確性で定義できる。化学で最も一般的に用いられているものは分子力学法に基づいており、粒子-粒子相互作用の古典的取扱い（構造変化やコンホメーション変化は再現できるが、化学反応は再現できない）を具体化している。

完全な量子力学的記述から古典的ポテンシャルへの簡略化は2つの主要な近似を伴う。1つ目はボルン＝オッペンハイマー近似である。この近似では電子のダイナミクスが非常に速く、核の運動に瞬間的反応すると考えることができる、と述べる。結果として、電子の動きと核の動きは別々に扱うことができる。2つ目の近似は、電子よりもかなり重い核を古典ニュートン動力学に従う点粒子として扱う。古典的分子動力学では、電子の影響は単一のポテンシャルエネルギー表面（通常は基底状態を表す）として近似される。

より細かい詳細が必要な時は、量子力学に基づくポテンシャルが用いられる。また、系の大部分を古典的に扱うが、化学的変換が起こる小さな領域を量子系として扱うハイブリッド古典/量子ポテンシャルも開発されている。

経験的ポテンシャル[編集]

化学で用いられる経験的ポテンシャルは力場と呼ばれることが多いのに対して、材料化学分野では原子間ポテンシャルと呼ばれる。

化学におけるほとんどの力場は経験的なものであり、化学結合と関連する結合力、結合角、結合二面角、ファンデルワールス力および静電価と関連する非結合力の和から成る。経験的ポテンシャルはアドホックな機能的近似によって限定的に量子力学的効果を表わす。これらのポテンシャルは原子電荷、原子半径の推定値を反映するファンデルワールスパラメータ、平衡結合長、結合角、結合二面角といった自由なパラメータを含む。これらは、詳細な電子構造（量子化学シミュレーション）あるいは弾性係数、格子パラメータ、分光測定といった経験的な物理的性質に対してフィッティングを行うことで得られる。

非結合性相互作用の非局所的な特性のため、これらは系の全ての粒子間の弱い相互作用を少なくとも含む。その計算は通常、MDシミュレーションの速度のボトルネックである。計算コストを下げるため、力場はシフト打ち切り半径、反応場アルゴリズム、粒子メッシュ・エバルト和、あるいはより新しい粒子-粒子-粒子-メッシュ（P3M）法といった数値的近似を用いる。

化学力場は一般にあらかじめ設定された結合様式を用いる（非経験的動力学法を除く）。したがって、化学力場は化学結合の切断の過程や反応を露にモデル化することができない。一方で、結合次数形式に基づいたもののような物理学におけるポテンシャルの多くは、系の複数の異なる接続や結合の切断を記述することができる^[11][12]。こういったポテンシャルの例としては、炭化水素のためのブレナー・ポテンシャル^[13]やそれをC-Si-H系^[14]とC-O-H系^[15]にさらに発展させたものがある。ReaxFFポテンシャル^[16]は、結合次数ポテンシャルと化学力場とを組み合わせた完全な反応力場と見なすことができる。

対ポテンシャルと多体ポテンシャル[編集]

非結合性エネルギーを表わすポテンシャル関数は、系の粒子間の相互作用全体の和として定式化される。多くの人気のある力場で採用されている最も単純な選択肢は、全ポテンシャルエネルギーが原子の対の間のエネルギー寄与の和から計算できる「対ポテンシャル」である。こういった対ポテンシャルの一例は非結合性レナード＝ジョーンズ・ポテンシャルであり、ファンデルワールス力を計算するために使われる。

${\displaystyle U(r)=4\varepsilon \left[\left({\frac {\sigma }{r}}\right)^{12}-\left({\frac {\sigma }{r}}\right)^{6}\right]}$

もう一つの例はイオン格子のボルン（イオン）モデルである。次式の第一項はイオンの対についてのクーロンの法則であり、第二項はパウリの排他原理によって説明される短距離反発であり、最終項は分散相互作用項である。大抵は、シミュレーションは双極子項のみを含むが、四極子項も同様に含まれることもある（バッキンガム・ポテンシャルとして知られる）。

${\displaystyle U_{ij}(r_{ij})={\frac {z_{i}z_{j}}{4\pi \epsilon _{0}}}{\frac {1}{r_{ij}}}+A_{l}\exp {\frac {-r_{ij}}{p_{l}}}+C_{l}r_{ij}^{-n_{l}}+\cdots }$

多体ポテンシャルにおいて、ポテンシャルエネルギーは互いに相互作用する3つ以上の粒子の効果を含む。対ポテンシャルを用いたシミュレーションでは、系の包括的な相互作用も存在するが、対ポテンシャル項を通じてのみ生じる。多体ポテンシャルにおいて、ポテンシャルエネルギーは原子の対全体の和によって表わすことができない。これは、これらの相互作用が高次項の組合せとして明確に計算されるためである。統計的見方では、変数間の依存性は一般に自由度の対ごとの積のみを用いて表現することはできない。例えば、炭素、ケイ素、ゲルマニウムのシミュレーションに元々使われ、その他の幅広い材料に対しても用いられているターソフ・ポテンシャル^[17]は3個の原子の群についての和を含む。このポテンシャルでは、原子間の角度が重要な要素である。その他の例としては、原子挿入法（EAM）^[18]や強結合二次モーメント近似（TBSMA）ポテンシャル^[19]がある。TBSMAポテンシャルでは、原子の領域における状態の電子密度は周囲の原子からの寄与の和から計算され、ポテンシャルエネルギー寄与はこの和の関数である。

半経験的ポテンシャル[編集]

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

強結合ポテンシャルとして知られる半経験的ポテンシャルには様々な種類があり、これらはモデル化される原子によって異なる。

分極可能なポテンシャル[編集]

ほとんどの古典的力場は分極率の効果を黙示的に含む（例えば量子化学計算から得られた部分電荷を拡大することによって）。これらの部分電荷は原子の質量に関して固定である。しかし、分子動力学シミュレーションはドルーデ粒子や変動電荷といった異なる手法を用いた誘導双極子の導入によって分極率を明示的にモデル化できる。これによって、局所的な化学的環境に応答する原子間の電荷の動的再分配が可能になる。

長年、分極可能MDシミュレーションは次世代シミュレーションとしてもてはやされてきた。水といった均一な液体については、分極率を含めることによって正確性の向上が達成されてきた^[20]。タンパク質についても有望な結果が得られている^[21]。しかしながら、シミュレーションにおいて分極率をどのように近似するのが最適化についてはいまだ不確かである^[要出典]。

ab-initio法におけるポテンシャル[編集]

古典的分子動力学では、単一のポテンシャルエネルギー表面（通常は基底状態）は力場によって表わされる。これはボルン＝オッペンハイマー近似の結果である。励起状態では、化学反応あるいはより正確な表現が必要な時は、電子の振る舞いを密度汎関数法といった量子力学的手法を用いることによって第一原理から得ることができる。これはab initio分子動力学（AIMD）と呼ばれる。電子の自由度を扱うコストから、このシミュレーションの計算コストは古典的分子動力学よりもかなり高い。これはAIMDがより小さな系あるいはより短い時間に制限されることを意味する。

Ab initio量子力学法は、トラジェクトリ中の配座について必要に応じてその場で系のポテンシャルエネルギーを計算するために使うことができる。この計算は反応座標の近傍で大抵行われる。様々な近似を使うことができるが、これらは経験的当て嵌めではなく理論的考察に基づいている。Ab-initio計算は、電子状態の密度やその他の電子的性質といった、経験的手法からは得ることのできない膨大な情報を与える。Ab-initio法を使用する大きな利点は、共有結合の切断あるいは形成を含む反応を調べる能力である。これらの現象は複数の電子状態に対応する。

ハイブリッドQM/MM法[編集]

QM（量子力学的）法は非常に強力である。しかしながら、その計算コストは高い。それに対してMM（古典的あるいは分子力学）法は高速だが、いくつかの制限がある（膨大なパラメータ化の必要性、得られたエネルギー推定値がそれ程正確ではないこと、共有結合が切断/形成する反応のシミュレーションに使うことができないこと、化学的反応に関する正確な詳細を与える能力に限界があること）。QM計算の利点（正確性）とMM計算の利点（速さ）を組み合わせた新たな手法が開発されている。これらの手法は混合あるいはハイブリッド量子力学/分子力学法（ハイブリッドQM/MM法）と呼ばれている^[22]。

ハイブリッドQM/MM法の最も重要な利点は速さである。最も分かりやすい場合において古典的分子動力学 (MM) を行うコストはO(n²) と見積られる（nは系中の原子の数）。これは主に静電相互作用項（全ての粒子がその他全ての粒子と相互作用する）のためである。しかしながら、打ち切り半径の使用、周期的対表の更新、粒子-メッシュ・エバルト (PME) 法の派生法によって、このコストをO(n) からO(n²) に減らすることができる。言い換えると、2倍の数の原子の系をシミュレーションすると、2倍から4倍の計算力を要することになる。一方で、最も単純なab-initio計算のコストは典型的にO(n³) あるいはそれ以上を見積られる（制限ハートリー＝フォック計算は～O(n^2.7) でスケールすることが示唆されている）。この制限を乗り越えるため、系の小さな部分が量子力学的に取り扱われ（典型的には酵素の活性部位）、残りの系が古典的に取り扱われる。

より洗練された実装では、QM/MM法は量子効果に対して敏感な軽い核（例えば水素）と電子状態の両方を扱うために存在する。これによって、水素の波動関数の生成を行うことができる。この方法論は、水素のトンネリングといった現象を調べるために有用である。QM/MM法が新たな発見をもたらした一つの例は、肝臓のアルコール脱水素酵素におけるヒドリド転移の計算である。この場合、水素原子のトンネリングが重要である（反応速度を決定する）^[23]。

粗視化表現[編集]

詳細なスケールの対極にあるのが、粗視化モデルと格子モデルである。系の全ての原子を露に（明示的に）表現する代わりに、ここでは原子の群を表現するために「擬原子」を用いる。非常に大きな系のMDシミュレーションは非常に大きな計算機資源を必要とするため、伝統的な全原子手法によって容易に調べることができない。同様に、長い時間スケール（1マイクロ秒を超える）の過程のシミュレーションは、多くの時間ステップを必要とするため、極めて計算コストが高い。これらの場合、粗視化（粗粒化）表現とも呼ばれる簡約表現を用いることによってこの問題に対処することができることもある。

粗視化（coarse graining、CG）手法の例としては、不連続分子動力学（CG-DMD）^[24][25]やGoモデルがある^[26]。粗視化はより大きな擬原子を用いることによって行なわれることもある。こういった合同原子（united atom）近似は生体膜のMDシミュレーションにおいて使用されてきた。電子的性質が興味の対象である系へのこういった手法の導入は、擬原子上の適切な電荷分布を使うことの困難さのため難しい^[27]。脂質の脂肪族末端は2から4のメチレン基を1つの擬原子としてまとめたいくつかの擬原子によって表わされる。

これらの非常に粗視的なモデルのパラメータ化は、モデルの挙動を適切な実験的データあるいは全原子シミュレーションへ合致させることによって、経験的に行われる。理想的には、これらのパラメータは自由エネルギーへのエンタルピー寄与とエントロピー寄与の両方を黙示的に考慮しなければならない。粗視化がより高い水準で行われる時、動力学的記述の正確性はより信頼できなくなる。しかし、よく粗視化されたモデルは、構造生物学、液晶の組織化、高分子ガラスの分野における幅広い疑問を調べるためにうまく使われてきている。

粗視化の応用の例を以下に挙げる。

• タンパク質のフォールディングの研究はアミノ酸毎に単一（あるいはいくつかの）擬原子を使ってしばしば行なわれる。
• 液晶の相転移は制限された幾何構造と異方性種を記述するGay-Berneポテンシャルを用いた計算の一方あるいは両方で調べられている。
• 変形中のポリマーガラスは、レナード＝ジョーンズポテンシャルによって記述され球を接続する単純な調和バネあるいは有限伸張性の非線形バネ (FENE; Finitely Extensible Nonlinear Elastic) を用いて研究されている。
• DNA超らせん化は塩基対当たり1-3の擬原子を用いて、またそれよりもさらに低い分解能で研究されている。
• 二重らせんDNAのバクテリオファージ内への詰め込みは二重らせんの1ターン（約10塩基対）を表わす1つの擬原子を使ったモデルによって調べられている。
• リボソームやその他の大きな系におけるRNA構造はヌクレオチド当たり1つの擬原子を用いてモデル化されている。

最も単純な粗視化の形は「合同原子（united atom）」であり、初期のタンパク質、脂質、核酸のMDシミュレーションのほとんどで使われた。例えば、CH₃メチル基の4原子全て（あるいはCH₂メチレン基の3原子全て）を露（明示的）に扱う代わりに、メチル基あるいはメチレン基全体を単一の擬原子によって表わす。この擬原子はもちろん、他の基とのファンデルワールス相互作用が適切な距離依存性を持つように適切にパラメータ化されなければならない。この種の合同原子の表現においては通常、水素結合に関与する能力のあるもの（極性水素）を除いて全ての明示的水素原子を消去する。この一つの例がCharmm 19力場である。

極性水素は通常モデルに保持される。これは水素結合の適切な取扱いが水素結合ドナー基とアクセプター基との間の指向性と静電相互作用のかなり正確な記述を必要とするためである。例えば水酸基は水素結合ドナーと水素結合アクセプターのどちらのなることができ、単一のOH擬原子ではこれを扱うことは不可能であろう。ここで留意すべきはタンパク質あるいは核酸中の原子の約半数は非極性水素であることであり、したがって合同原子を使用することによって計算時間を相当短縮することができる。

操舵分子動力学 (SMD)[編集]

操舵分子動力学（Steered molecular dynamics; SMD）シミュレーションでは、望む自由度に沿ってタンパク質の構造を引っ張ることによってその構造を操作するためにタンパク質に力を印加する。これらの実験は原子レベルでのタンパク質における構造変化を明らかにするために用いることができる。SMDは機械的な折り畳み構造のほどけや伸長といった出来事をシミュレーションするためにしばしば用いられている^[28]。

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

応用例[編集]

分子動力学は多くの科学分野で使われている。

• 単純化された生物学的折り畳み過程の最初のMDシミュレーションは1975年に発表された。Nature誌で発表されたそのシミュレーションは現代のタンパク質折り畳み計算の広大な領域への道を開いた^[32]。
• 生物学的過程の最初のMDシミュレーションは1976年に発表された。Nature誌で発表されたそのシミュレーションはタンパク質の運動が単なる飾りではなく機能に必須であることの理解への道を開いた^[33]。
• MDはheat spike regimeにおける衝突カスケード、すなわちエネルギー中性子とイオン放射が固体および固体表面上で持つ効果を取り扱うための標準的手法である^[34][35]。
• MDシミュレーションはゴーシェ病の原因である最も一般的なタンパク質変異N370Sの分子基盤を予測することにうまく応用された^[36]。後続の論文では、これらの目隠し予測が同じ変異についての実験結果と驚く程に高い相関を見せることが示された^[37]。
• MDシミュレーションは金属表面上の水薄膜の分離圧に対する表面電荷の影響について調べるために用いられている^[38]。
• MDシミュレーションは透過型電子顕微鏡の画像特徴を理解するためにマルチスライス画像シミュレーションと共に用いられる^[39]。

以下の生物物理学的例は非常に大きな系（完全なウイルス）あるいは非常に長いシミュレーション時間（1.112ミリ秒まで）のシミュレーションを行うための注目に値する成果を示している。

• 完全なサテライトタバコモザイクウイルス（STMV）のMDシミュレーション（2006年、規模: 100万原子、シミュレーション時間: 50ナノ秒、プログラム: NAMD）。このウイルスは小さい20面体植物ウイルスであり、タバコモザイクウイルス (TMV) による感染の症状を悪化させる。分子動力学シミュレーションは、ウイルス集合の機構を詳細に調べるために用いられた。全STMV粒子はウイルスカプシド（被覆）を作り上げる単一タンパク質の同一の複製物60個と1063ヌクレオチドの一本鎖RNAゲノムから構成される。1つの重要な発見は、RNAが内部にない時はカプシドが非常に不安定であるということである。このシミュレーションは2006年のデスクトップコンピュータ1台では完了するのに約35年を要する。したがって、シミュレーションは並列に接続した多数のプロセッサによって行われた^[40]。
• ビリンタンパク質の頭部断片の全原子力場による折り畳みシミュレーション（2006年、規模: 2万原子、シミュレーション時間: 500マイクロ秒、プログラム: Folding@home）。このシミュレーションは参加した世界中のパーソナルコンピュータの20万CPU上で実行された。これらのコンピュータにはFolding@homeプログラムがインストールされていた。ビリン頭部タンパク質の動力学的特性は、連続したリアルタイムコミュニケーションを行わないCPUによる多くの独立した短時間のシミュレーションを用いることによって詳細に調べられた。使われた1つの手法が、特定の開始コンホメーションの折り畳みがほどける前の折り畳みの確率を測定するPfold値解析である。Pfoldは遷移状態構造と折り畳み経路に沿ったコンホメーションの規則化に関する情報を与える。Pfold計算におけるそれぞれのトラクジェクトリは比較的短くてもよいが、多くの独立したトラクジェクトリが必要である^[41]。
• 長い連続トラクジェクトリシミュレーションが、超並列スーパーコンピュータアントン上で実行された。発表された最長のアントンを用いて実行されたシミュレーション結果は355 KにおけるNTL9の1.112ミリ秒シミュレーションである。2番目は、同じ構造について独立して行われた1.073ミリ秒シミュレーションである^[42]。『How Fast-Folding Proteins Fold』において、研究者のKresten Lindorff-Larsen、Stefano Piana、Ron O. Dror、David E. Shawは「12種類の構造的に多様なタンパク質の折り畳みに内在する一連の一般原理を明らかにする100 μ秒から1 m秒の間の範囲に渡る原子レベルでの分子動力学シミュレーションの結果」について議論した。専用のカスタムハードウェアによって可能になったこれらの多様な長いトラクジェクトリの調査から、彼らは「ほとんどの場合において、折り畳みは、非折り畳み状態において形成される要素の傾向と高度に相関した順序でネイティブ構造の要素が現われる単一の支配的経路を取る」と結論付けた^[42]。別の研究において、300 Kにおけるウシ膵臓トリプシンインヒビター（BPTI）のネイティブ状態動力学の1.031ミリ秒シミュレーションを行うためにアントンが使われた^[43]。
• これらの分子シミュレーションは、材料除去の機構や道具の形状、温度、切断速度や切断力といった加工パラメータの影響について理解するために用いられている^[44]。また、数層のグラフェン^[45][46]やカーボンナノスクロールの剥離の背後にある機構を調べるためにも用いられた。

分子動力学アルゴリズム[編集]

積分器[編集]

• シンプレクティック積分器
• ベルレ＝ストーマー積分
• ルンゲ＝クッタ積分
• ビーマンのアルゴリズム
• 拘束アルゴリズム（拘束された系）

短距離相互作用アルゴリズム[編集]

• セル・リスト
• ベルレ・リスト
• 結合相互作用

長距離相互作用アルゴリズム[編集]

• エバルト和
• 粒子メッシュエバルト（Particle mesh Ewald、PME）
• 粒子–粒子–粒子–メッシュ（P3M）
• Shifted force method
• Isotropic periodic sum method

並列化戦略[編集]

• 領域分割法

分子動力学シミュレーションソフトウエアパッケージ[編集]

• Advance/PHASE
• AMBER
• CHARMM
• Desmond
• DL_POLY
• GROMACS
• LAMMPS
• Materials Studio
• NAMD
• PEACH
• PRESTO
• Quantum ESPRESSO
• SCIGRESS
• VASP
• Winmostar

脚注[編集]

参考文献[編集]

関連項目[編集]

• 2体ポテンシャル
• 3体ポテンシャル
• 多体ポテンシャル
• 多体問題
• セル多重極子展開法
• ベルレ法
• ギア法
• 原子挿入法 (EAM)
• 量子分子動力学
• 反対称化分子動力学
• 計算物理学
• 物性物理学
• 計算化学