MD計算における温度制御について

MD計算の目的の一つに指定された熱力学条件における系の性質を解析・予測することが挙げられます。特定の温度での状態をシミュレーション上で実現するために，系の自由度に対する熱浴からの影響を何かしらの形でモデル化します。そのモデル化には様々な手法がこれまで提案されており，MD計算エンジンによっては複数手法が実装されていてどれを使うべきか迷ってしまう方もいるかと思います。この記事では，幾つかの温度制御法の特徴と注意点を解説したいと思います。

Global thermostatとLocal thermostat

温度制御法は大きく分けて，Global thermostatとLocal thermostatがあります。
Global thermostatは系内の全運動エネルギーを制御することで対象温度での状態を実現します。代表的なGlobal thermostatとしてNosé–Hoover法やBussi法などが挙げられます。
一方，Local thermostatは系内の個々の運動エネルギーを制御することで対象温度での状態を実現します。代表的なLocal thermostatとしてはLangevin法などが挙げられます。
一見，理論的にはGlobal thermostatとLocal thermostatのどちらも使用しても問題ないように思えるのですが，実用上はどちらも懸念点があるため，どのthermstatを選択するかは重要となる場合があります。

Global thermostatの問題点

Global thermostatは全運動エネルギーをモニタリングして対象温度に制御するため，個々の運動エネルギーの分布に対しては特別な制御がなされません。その結果，異なる自由度が同一の温度に制御されないという事態が生じえます。例えば，トータルだと300Kになっているが，タンパク温度が290Kで溶媒温度が305Kになっているような状況です。世間的にはHot-Solvent/Cold-Solute Problemと称されたりします。 このような状態が瞬間的ではなくシミュレーション中に持続することになると，溶媒⇒タンパクのheat flowが恒常的に存在することになりますので，厳密には平衡状態と言えないことになります。また，300Kでのタンパクの状態を調べたいのに，シミュレーション上では290Kの状態が実現されているというミスマッチも生じていることになります。
「異なる自由度が同一の温度に制御されない」という問題は，究極的には系全体の並進/回転に全運動エネルギーが集中して内部自由度が凍結してしまう事態にも発展しえます。これは，Flying ice cube effectと称されます。WikiではBussi法ではこの問題が生じないと断言されていますが，少なくとも自分の経験ではNosé–Hoover法やBussi法のどちらでもFlying ice cube effectが生じたことがあります。

Local thermostatの問題点

Local thermostatは個々の運動エネルギーの分布に対して温度制御がなされるため，異なる自由度が同一の温度に制御されないという事態は生じません。その結果，Global thermostatで生じうるHot-Solvent/Cold-Solute ProblemやFlying ice cube effectの心配はありません。ですが，特に樹脂や生体分子といったポリマー分子を扱うMD計算では，時間スケールの遅いモーションが過剰に抑圧されてしまう問題が散見されます。これは，個々の原子に温度制御のための力が課されるために，集団的な動きをするモードが生じにくい状態になっていることに由来します。その結果，Hot-Solvent/Cold-Solute Problemが生じているGlobal thermostatでのサンプリングよりも，サンプリング効率が著しく低下してしまっているという問題が生じえます。

Multiple global thermostat

Global thermostatとLocal thermostatはどちらも問題点があるのですが，生体分子系のMD計算ではLocal thermostatの問題の方が深刻です。また，Hot-Solvent/Cold-Solute Problemもできれば避けたい問題です。そのため，生体分子系のMD計算では，サンプリングのMD計算時には幾つかの部分構造（タンパク，溶媒，etc）にそれぞれGlobal thermostatによる温度制御を実施することが多い印象です。これは，Global thermostatとLocal thermostatの中間をとっているような状態です。ですが，たった２つのthermostatでもlocal thermostatと同様な問題が生じることが指摘されており，この手法が最もベストとも言えないことには留意すべきかと思います。

NVE条件下でのサンプリング

Global thermostatとLocal thermostatどちらの問題も，起源は「熱浴からの影響を何かしらの形でモデル化」したことに由来します。そのため，そのモデル化の要素を排除したNVE条件下では上記のいずれの問題も生じません。NVE条件下でのMD計算をNVTやNPT条件下でのサンプリングと結びつける方法は↓の記事をご参照ください。

NVE条件下でのサンプリングのデメリットは，計算コストがより多くかかってしまう点です。NVEの「E」に関してNVTやNPTより厳しい計算精度が要求されますし，複数の「E」に対してMD計算が必要な場合も考えられます。

結局どうすればよいのか？

状況に応じて使い分けるというのが良いと思います。
計算対象がフラグメントサイズの低分子＋水溶媒からなる系でしたら，Langevinの適用は問題にならないと思います。
生体分子シミュレーションのようにタンパク＋水溶媒（必要に応じて＋脂質膜）から構成され複数の時間スケールの異なる自由度が混在する系の場合，できる限りthermostatによるアーチファクトをなくしたいのであれば，計算コストをかけてでもNVEでサンプリングすべきです。ある程度のアーチファクトを覚悟の上で（multiple）Global thermostatを適用するのも考え方の一つとしてあると思います。少なくとも樹脂や生体分子といったポリマー分子を扱うMD計算では，サンプリングの際にLangevinを使用するのは推奨しません。