ASEでFixAtomsを用いて、一部の原子を固定してNPTclassを実行してNVTを実行した際、設定した温度より計算される温度が高くなる現象が報告されています。この現象の原因と対応策を紹介します。
本記事のコードは、ASE からの抜粋です。コードの著作権は © Copyright 2025, ASE-developers に帰属し、GNU Lesser General Public License version2.1 の条項に基づいて公開されています。
■原因
分子動力学(MD)シミュレーションの温度は各原子の持つ運動エネルギーにより決定されます。NPTclassでは、最初に、到達させたい運動エネルギーの総量を、固定原子を除外せずに計算しています。
これはFixAtomsで固定された原子も周囲と同じ運動量を持っていることを仮定していますが、固定された原子がある場合、の運動エネルギーはゼロとして計算するのが正しく、この仮定が成り立ちません。
そのため、所望の温度にならない現状が生じております。
ソースコードを見ながら、詳細を以下で説明します。
NPTclassの_calculateconstantsで設定した温度の運動エネルギーは以下のように計算されています。
self.desiredEkin = 1.5 * (n - 1) * self.temperaturenは原子数です。系の自由度は3nですが、系全体の並進がないために自由度が-3されます。そのため、ある温度下の運動エネルギーは0.5 * (3n-3) * temperatureとなり、上の実装と一致します。自由度が原子数を指定しているため、これはFixAtomsの効果を反映していません。
そしてNPTclassのrunでは以下のように未来の位置q_futureと過去の位置q_pastから運動量を計算しています。
self.atoms.set_momenta(
np.dot(self.q_future - self.q_past, self.h / (2 * dt)) *
self._getmasses())
このとき、FixAtomsで原子が固定されていると、その原子の運動量は0になり、Atomsのget_get_kinetic_energyではFixAtomsの効果を反映した運動エネルギーを計算することになります。
そうすると、設定した温度から期待される系の運動エネルギー(self.desiredEkin)とFixAtomsの効果を加味した実際の運動エネルギー(self.atoms.get_kinetic_energy())に差が生じ、deltazetaが負になり、運動エネルギーを大きくしようとします。
deltazeta = 2 * dt * self.tfact * (self.atoms.get_kinetic_energy() -
self.desiredEkin)しかし、Atomsのget_temperatureでは、系の自由度を加味して温度を計算する実装になっています。そのため、目標の運動エネルギーが、入力したtemperatureで想定されている値よりも、高く設定されているので、系の実際の温度も高くなっています。表示される温度は系の温度を反映しています。
def get_temperature(self):
"""Get the temperature in Kelvin."""
ekin = self.get_kinetic_energy()
return 2 * ekin / (self.get_number_of_degrees_of_freedom() * units.kB)■対応策
ここでは、3つの対応策を紹介します。
①設定温度を補正してNPTclassを使用する
以下のように設定温度を手動で補正していただくことで設定した温度に近づけることができます。 補正温度 = 設定温度*(3*[自由な原子数] - 3) / (3*[自由な原子数] + 3*[固定した原子数] - 3)
②Langevinclassに変更する
LangevinclassはFixAtomsと一緒にお使いいただけます。
③NVTBerendsenclassに変更する
NVTberendsenclassもFixAtomsに対応しているため、ご利用いただけます。
ただし、正しい熱平衡状態を生成する保証はありませんので、別の温度制御法のご利用をお勧めします。