はじめに¶

計算化学と分子モデリング¶

GROMACS は、分子動力学シミュレーションとエネルギー最小化を実行するためのエンジンです。これらは、計算化学および分子モデリングの分野における多くの技術の1つです。*計算化学*とは、分子の量子力学から、大規模で複雑な分子集団のダイナミクスまで、化学における計算技術の使用を示す名称です。*分子モデリング*とは、現実的な原子モデルを用いて複雑な化学システムを記述する一般的なプロセスであり、詳細な原子レベルでの知識に基づいて、マクロな特性を理解し予測することを目的としています。分子モデリングは、特に、現実的なシステムの正確な物理特性を予測する必要がある新しい材料の設計によく使用されます。

マクロな物理的特性は、以下のように区別できます。

静的平衡の特性、例えば、酵素に対する阻害剤の結合定数、系の平均ポテンシャルエネルギー、または液体のラジアル分布関数など。
動的な特性または非平衡状態の特性、例えば液体の粘度、膜中の拡散プロセス、相変化のダイナミクス、反応速度、または結晶中の欠陥のダイナミクスなど。

選択する手法は、質問の内容と、現在の最先端技術において、その手法が信頼できる結果をもたらす可能性によって決まります。理想的には、（相対論的な）時間依存のシュレーディンガー方程式は、分子系の特性を高い精度で記述しますが、数個の原子の平衡状態よりも複雑なものは、この*ab initio*レベルでは扱えません。したがって、近似が必要となります。システムの複雑さと、関心のあるプロセスの時間範囲が長いほど、必要な近似はより厳密になります。ある時点で（通常よりも大幅に早い段階で）、*ab initio*アプローチは、使用しているモデルの*経験的*パラメータ化によって補強または置き換える必要があります。原子相互作用の物理学に基づいたシミュレーションが、システムの複雑さのために失敗する場合、分子モデリングは、既知の構造および化学データとの類似性分析に基づいています。QSAR（定性構造活性関係）法や、多くのホモロジーに基づくタンパク質構造予測は、このカテゴリーに属します。

マクロな特性は、常に代表的な統計的な系（平衡状態または非平衡状態）の分子系の集合平均として表されます。分子モデリングにおいては、このことは以下の2つの重要な意味を持ちます。

単一の構造、たとえそれがグローバルなエネルギー最小構造であっても、その知識だけでは不十分です。マクロな特性を計算するためには、特定の温度における代表的な構造の集合を生成する必要があります。しかし、これは、自由エネルギーに基づく熱力学的な平衡特性、例えば相平衡、結合定数、溶解度、分子構造の相対的な安定性などを計算するために不十分です。自由エネルギーと熱力学ポテンシャルの計算には、分子シミュレーション技術の特別な拡張が必要です。
分子シミュレーションは、原理的には構造と運動の詳細を提供しますが、これらの詳細は、一般的に興味のあるマクロな特性には関連性がありません。これにより、相互作用の説明を簡略化し、無関係な詳細を平均化することができます。統計力学の科学は、このような簡略化のための理論的な枠組みを提供します。原子を1つの単位として扱うことから、集団座標の少ない数で運動を記述し、平均力ポテンシャルと確率的ダイナミクスを組み合わせた溶媒分子の平均化、さらにはref:9。

代表的な平衡状態のアンサンブルを生成するために、以下の2つの方法が利用可能です:

*モンテカルロシミュレーション*と
分子動力学シミュレーション

非平衡エンsemblesの生成やダイナミックな現象の解析には、2番目の方法のみが適しています。モンテカルロシミュレーションは、MD（力の計算が不要）よりも単純ですが、同じ計算時間でMDと同程度の統計的精度を得ることはできません。したがって、MDの方が汎用的な手法です。初期配置が非常に非平衡の場合、力は過度に大きくなり、MDシミュレーションが失敗する可能性があります。そのような場合は、堅牢な*エネルギー最小化*が必要です。エネルギー最小化を行うもう一つの理由は、システムの全運動エネルギーを取り除くことです。複数の「スナップショット」をダイナミックシミュレーションから比較する必要がある場合、エネルギー最小化は構造とポテンシャルエネルギーの熱ノイズを低減し、それらをより適切に比較することができます。

分子動力学シミュレーション¶

MDシミュレーションは、相互作用する \(N\) 個の原子のニュートンの運動方程式を解きます。

(1)¶\[m_i \frac{\partial^2 \mathbf{r}_i}{\partial t^2} = \mathbf{F}_i, \;i=1 \ldots N.\]

これらの力は、ポテンシャル関数 :math:`V(mathbf{r}_1, mathbf{r}_2, ldots, mathbf{r}_N)`の負の導関数です。

(2)¶\[\mathbf{F}_i = - \frac{\partial V}{\partial \mathbf{r}_i}\]

これらの方程式は、小さな時間ステップで同時に解かれます。システムの状態をある程度の時間追跡し、温度と圧力が適切な値に維持されるように制御し、座標を定期的に出力ファイルに書き込みます。時間に対する座標は、システムの*軌跡*を表します。初期の変化の後、システムは通常、*平衡状態*に達します。平衡状態の軌跡を平均化することで、出力ファイルから多くのマクロな特性を抽出できます。

この時点で、MDシミュレーションの限界を考慮することが有益です。ユーザーはこれらの限界を認識し、シミュレーションの精度を評価するために、既知の実験的特性について常に確認を行う必要があります。以下に近似をリストします。

シミュレーションは古典的なもの

ニュートンの運動方程式を使用することは、原子の運動を記述するために*古典力学*を使用することを意味します。これは、通常の温度でほとんどの原子にとっては問題ありませんが、例外もあります。水素原子は非常に軽いため、プロトンの運動は時々重要な量子力学的な性質を示します。たとえば、プロトンは水素結合を通過する際に、ポテンシャル障壁を*トンネル*することができます。このようなプロセスは、古典的なダイナミクスでは適切に扱えません！低温のヘリウム液体は、古典力学が破綻する別の例です。ヘリウム自体は私たちにとって直接的な問題ではないかもしれませんが、コバレント結合の高い周波数振動は私たちを心配させるべきです！古典的な調和振動子の統計力学は、共鳴周波数 \(\ u\) が \(k_BT/h\) に近づいたり超えたりした場合、実際の量子振動子のものとは顕著に異なります。室温では、周波数 \(\sigma = 1/\lambda = \ u/c\) で \(h \ u = k_BT\) が約 200 cm\(^{-1}\) になります。したがって、たとえば 100 cm\(^{-1}\) よりも高いすべての周波数は、古典的なシミュレーションで問題を引き起こす可能性があります。これは、ほとんどの結合と結合角の振動が疑わしく、水素結合による翻訳または回転的なH-結合振動も古典的な限界を超えていることを意味します（参照：Table 1）。私たちは何をすべきでしょうか？

表 1 典型的な振動周波数（波数）は、分子および水素結合性液体に見られる。比較：`kT/h = 200~mathrm{cm}^{-1}` (300 K)。¶
債のタイプ	振動の種類	波数 \(\mathrm{cm}~^{-1}\)
C-H、O-H、N-H	伸張	3000～3500
C=C, C=O	伸張	1700--2000
HOH	曲がり	1600
C-C	伸張	1400～1600
H\(_2\)CX	スイス、ロック	1000～1500
CCC	曲がり	800--1000
O-H\(\cdots\)O	振動	400～700
O-H\(\cdots\)O	伸張	50～200

さて、実質的な量子・動的シミュレーション以外にも、以下の2つの方法があります。
1. もし、ハarmonicなオシレーターを使用して分子動力学シミュレーションを行う場合、以下の修正を行う必要があります：総内部エネルギー：U = E_{kin} + E_{pot} および比熱：C_V`（およびエンタルピー：`S および自由エネルギー：A または G が計算されている場合） 1次元オシレーター（周波数：u）のエネルギーおよび比熱に対する修正は、以下の通りです：<refMcQuarrie76>
  
  (3)¶\[U^{QM} = U^{cl} +kT \left( {\frac{1}{2}}x - 1 + \frac{x}{e^x-1} \right)\]
  
  (4)¶\[C_V^{QM} = C_V^{cl} + k \left( \frac{x^2e^x}{(e^x-1)^2} - 1 \right)\]
  
  ここで、\(x=h\ u /kT\)。古典的な振動子は、過剰なエネルギー (\(kT\)) を吸収しますが、高周波の量子振動子は、ゼロ点エネルギーレベルである \(\frac{1}{2} h\ u\) で基底状態にあります。
2. 方程式における結合（および結合角）を*制約*として扱うことができます。その理由は、量子オシレーターが基底状態にある場合、古典的なオシレーターよりも結合の制約に近い振る舞いをすることです。この選択の合理的な実用的な理由は、最も高い周波数が除去された場合に、アルゴリズムがより大きな時間ステップを使用できることです。実際には、結合を制約する場合とオシレーターとして扱う場合で、時間ステップを4倍にすることができます：<refGunsteren77>。|Gromacs|は、結合と結合角に対してこのオプションを提供しています。後者の柔軟性は、現実的な運動と構成空間のカバーを可能にするために非常に重要です：<refGunsteren82>。

電子は基底状態にある: MDシミュレーションでは、原子の位置のみに基づいて定義された*保守的な*力場を使用します。これは、電子の動きを考慮しないことを意味します。電子は、原子の位置が変化したときに、瞬時にそのダイナミクスを調整し（ボルン・オペニハイマー近似）、基底状態に留まることが想定されています。これは、ほとんどの場合、問題ありません。ただし、電子移動や電子的に励起された状態は、適切に扱ることができません。また、化学反応も適切に扱ることができませんが、現時点では反応を避ける理由がいくつかあります。
力場は近似的なものである: 力場は力を提供します。それらはシミュレーション手法の一部ではなく、必要に応じてまたは知識が向上するにつれて、ユーザーによってパラメータが変更できます。ただし、特定のプログラムで使用できる力の形式には制限があります。|Gromacs|に組み込まれている力場は、第4章で説明されています。現在のバージョンでは、力場はペアに依存したものであり（長距離のクーロン力を除く）、極性度を考慮することができず、結合相互作用の微調整は含まれていません。そのため、以下のリストにはいくつかの制限を含める必要があります。それ以外は、生物学的に関連する分子を水溶液で使用する場合に、非常に有用で信頼性が高いです！
力場はペア依存: これは、すべての*非結合*力は、非結合ペア相互作用の合計から生じることを意味します。最も重要な例である原子の極性による相互作用のような、ペアに依存しない相互作用は、*有効なペアポテンシャル*で表現されます。平均的な非ペア依存的な寄与のみが組み込まれます。これは、ペア相互作用が純粋ではない、つまり、テストシステムでパラメータ化されたモデルとは異なり、孤立したペアや、それと著しく異なる状況には適用されないことを意味します。実際には、有効なペアポテンシャルは実用上問題ありません。しかし、極性による寄与を無視することで、原子中の電子が本来持つ誘電定数を提供しません。たとえば、実際の液体アルカンは誘電定数がわずか2を超えており、これにより(部分的な)電荷間の長距離静電相互作用が減少します。したがって、シミュレーションは長距離のクーロン項を過大評価します。幸いなことに、次の項目がこの効果をある程度補います。
長距離の相互作用は中断されます: GROMACS では、通常、レンナー・ジョーンズ相互作用、および場合によってはクーロン相互作用にもカットオフ半径を使用します。「最小画像法」を使用する GROMACS では、周期境界条件における各粒子のうち、1つの画像のみがペア相互作用に使用されるため、カットオフ半径はボックスサイズの半分を超えてはなりません。これは、大規模なシステムでは依然として大きいため、特に電荷を持つ粒子を含むシステムでは問題が発生する可能性があります。しかし、電荷がカットオフ境界に集中したり、非常に不正確なエネルギーが得られたりするなど、深刻な問題が発生する可能性があります。このようなシステムでは、粒子メッシュ・エワルドなどの実装されている長距離静電相互作用アルゴリズムを使用することを検討してください：<refDarden93>、<refEssmann95>。
境界条件は不自然: システムサイズが小さい（例えば、10万個の粒子でも小さい）場合、粒子群は周囲（真空）との間で多くの不要な境界を持つことになります。このような状態を避け、実物の相境界をシミュレーションしたい場合は、周期境界条件を使用する必要があります。液体は結晶ではないため、不自然なものが残ります。この項目は最後に記述していますが、これは最も良い解決策の一つです。大規模なシステムでは誤差は小さいですが、内部の空間的な相関が強い小さなシステムでは、周期境界条件が内部の相関を増強する可能性があります。その場合は、システムサイズの影響を注意深く観察し、テストする必要があります。特に、長距離静電力を計算する際に、格子和を使用する場合は、この点に注意が必要です。なぜなら、格子和は、場合によっては不必要な秩序を導入することが知られているからです。

エネルギー最小化と探索手法¶

上記のように、:ref:`Compchem`で述べたように、多くの場合はエネルギー最小化が必要です。 |Gromacs|は、:ref:`EM`で詳細に説明されているように、ローカルエネルギー最小化のためのいくつかの方法を提供します。

(マクロ)分子系のポテンシャルエネルギー関数は、非常に複雑な地形（または*ハイパーサーフェス*）であり、多数の次元にわたるものです。存在するのは最も深い点、すなわち*グローバル・ミニマム*、そして非常に多くの*局所・ミニマム*です。これらの局所・ミニマムでは、ポテンシャルエネルギー関数の座標に対するすべての微分がゼロになり、すべての二階微分が非負になります。二階微分をまとめた行列、すなわち*ヘッセ行列*は、非負の固有値を持ちます。固有値がゼロになるのは、翻訳と回転に対応する、分子の集合座標のみです。局所・ミニマムの間には、ヘッセ行列に一つの負の固有値しか持たない*鞍点*が存在します。これらの点は、システムが一つの局所・ミニマムから別の局所・ミニマムへ移動するための山越え地点です。

すべての局所的な最小値、包括的な最小値を含む、およびすべての鞍点を知ることで、関連する構造と配置、およびそれらの自由エネルギー、ならびに構造変化のダイナミクスを記述できるようになります。しかし、構成空間の次元と局所的な最小値の数は非常に高いため、完全な調査を得るために十分な数の点から空間をサンプリングすることは不可能です。特に、実用的な時間内でグローバルな最小値を保証するような最適化方法は存在しません。実用的な方法が存在しますが、その速度は様々です：<refGeman84>。ただし、初期配置が与えられた場合、最も近い局所的な最小値を求めることができます。「最も近い」という言葉は、この文脈では必ずしも「幾何学的に最も近い」という意味ではありません（つまり、座標の差の最小和ではありません）。これは、最も急な局所的な勾配を下りながら、体系的に移動することで到達できる最小値を意味します。この最も近い局所的な最小値を特定することが、|Gromacs|があなたのためにできることのすべてです。グローバルな最小値を特定し、その過程で他の最小値を発見したい場合は、温度と結合されたMDを試すことが最善の方法です。システムをある程度の時間、高い温度で実行し、その後、必要な温度にゆっくりと冷却します。これを繰り返してください！もし、融点またはガラス転移温度が存在する場合は、その温度よりもわずかに低い温度でしばらく保持し、巧妙な方法でゆっくりと冷却することが賢明です。これは「シミュレーションアンニリング」と呼ばれます。物理的な制約は必要ありませんので、このプロセスを加速するために想像力を活用できます。よく機能するトリックの1つは、水素原子をより重くすることです（質量は10など）：これにより、水素原子の通常非常に速い動きが遅くなることはありませんが、システム内のより遅い動きにはほとんど影響を与えません。また、時間ステップを3倍または4倍にすることができます。また、検索手順中にポテンシャルエネルギー関数を修正することもできます。たとえば、障壁を取り除く（二面角関数を削除するか、反発ポテンシャルをソフトコアポテンシャルに置き換える<refNilges88>）などです。ただし、常に正しい関数をゆっくりと復元するように注意してください。最も大きな構造変化を可能にする最適な検索方法は、4次元空間への探索を許可すること<refSchaik93>です。ただし、これは|Gromacs|の標準機能を超えた追加のプログラミングが必要です。

以下の3つのエネルギー最小化手法があります。

これらは、単に関数評価のみを必要とするものです。例えば、シンプソックス法とそのバリエーションがあります。ステップは、これまでの評価の結果に基づいて決定されます。もし微分情報が利用できる場合は、この情報を利用する手法よりも劣ります。
これらの方法は、ポテンシャルエネルギーのすべての座標に対する偏微分がMDプログラムで既知である（これは力の負の値に等しい）ため、MDプログラムの修正に適しています。
これらは、2次微分情報を使用する方法も含まれます。これらの方法は、最小点付近での収束特性に優れています。1回のステップで二次ポテンシャル関数を最小化できます！ただし、:math:`N`個の粒子に対して、:math:`3Ntimes 3N`の行列を計算、保存、および逆行列化する必要があります。2次微分を求めるための追加のプログラミングに加えて、ほとんどの興味のあるシステムでは、これには十分な計算能力がありません。中間的な方法として、ハッセ行列を動的に構築する方法もありますが、これらも過剰なストレージ要件に悩まされます。そのため、|Gromacs|はこれらの方法を避けます。

GROMACS で利用できる 急降下法 は、2 階のアルゴリズムです。これは、前のステップで蓄積された履歴を考慮することなく、単に負の勾配（つまり、力の方向）の方向にステップを進めます。ステップサイズは調整され、探索が高速になるようにしますが、常に下り方向に移動します。これは、シンプルで堅実な方法ですが、やや非効率的です。特に、局所的な最小値の近くでは、収束が非常に遅くなることがあります！収束が速い 共役勾配法 (例: 19) は、前のステップからの勾配情報を利用します。一般的に、急降下法は、局所的な最小値に非常に速く近づくことができますが、共役勾配法は、局所的な最小値に 非常に 近づくことができますが、最小値から離れた場所では性能が低下します。GROMACS は、L-BFGS 最適化アルゴリズムもサポートしており、共役勾配法とほぼ同等の性能を発揮しますが、特定の状況ではより高速に収束することがあります。