樊哲勇《分子动力学模拟》python实例 | 4.用简谐振子展示控温算法的效果

樊哲勇老师最近出版的书籍《分子动力学模拟》中，使用C++语言进行分子模拟程序的实现。未了更好地理解分子动力学模拟，AI4S1O1-材料小队学习志愿者决定针对樊老师书籍中的内容，进行python语言的程序实现以及注释说明。

本节内容对应樊老师书中的第四章

樊哲勇《分子动力学模拟》python实例 | 概述

本章目标

1. 了解分子动力学模拟中的温度控制常用算法
2. 进行python编码实现
3. 进行不同算法之间的对比
4. 控温知识扩展

一、控温算法

1.1 定义

分子动力学（MD）控温算法（thermostat）是指在MD模拟中用来控制系统温度的算法。由于MD模拟是通过数值积分牛顿运动方程来计算粒子的运动轨迹，这种方法默认的是微正则系综（NVE系综），即系统的能量和粒子数目是固定的。在实际研究中，我们通常希望在正则系综（NVT系综，温度固定）或等压等温系综（NPT系综，温度和压力固定）下进行模拟，因此需要引入控温算法来调节系统温度。

使用控温算法的原因：

1. 研究温度依赖性现象：许多物理和化学现象都依赖于温度，例如相变、反应速率、扩散等。在研究这些现象时，需要精确控制系统温度。

2. 模拟实际实验条件：实验室中的大多数实验都是在恒温条件下进行的，通过控温算法，可以更好地与实验结果进行对比和验证。

3. 系统稳定性：在长时间模拟中，系统的温度可能会因为能量的积累或损失而发生变化，控温算法可以帮助维持系统在目标温度附近，确保模拟结果的可靠性。

1.2 常见的控温算法

Berendsen Thermostat： 通过将系统的温度逐渐调整到目标温度，调节系统的温度变化速率。优点是简单且收敛速度快，缺点是不能准确生成正则系综。

Nose-Hoover Thermostat：通过引入虚拟热浴和额外自由度来控制系统温度，可以准确产生正则系综。缺点是可能在小系统或短时间内产生温度振荡。

Langevin Thermostat：通过引入随机力和阻尼力来模拟与外界环境的热交换。优点是能够有效控制温度，并且在处理溶液和生物大分子时非常有用。

Andersen Thermostat：通过随机选择粒子进行速度重置，从而控制温度。优点是能够很好地热化系统，但对动力学过程的影响较大。

二、控温算法代码实现

Berendsen 控温算法

Berendsen 控温算法的核心思想是通过一个比例反馈机制来调整系统的动能，使系统的温度逐渐接近目标温度。具体来说，系统的温度 𝑇和目标温度 𝑇0的关系可以通过以下公式进行调整： $$\frac{dT}{dt}=\frac{T_0-T}{\tau }$$ 其中，𝜏是一个时间常数，决定了温度调节的速度。通过这个公式，我们可以得到新的温度 𝑇′： $$T^{\prime }=T+\frac{T_0-T}{\tau }\Delta t$$ 为了实现温度的调节，我们需要调整系统中每个粒子的速度 𝑣𝑖。调整后的速度vi′可以通过以下公式计算： $$v_i^{\prime }=\lambda v_i$$ 其中，𝜆是一个缩放因子，用于调整粒子的速度： $$\lambda =\sqrt{1+\frac{\Delta t}{\tau }\left(\frac{T_0}{T}-1\right)}$$ 这里，Δ𝑡是时间步长，𝑇是当前系统的温度，𝑇0是目标温度。

可知算法步骤：

1. 计算当前系统温度 𝑇：通过系统的动能计算当前温度。
2. 计算缩放因子 𝜆：使用上面的公式计算速度调整因子。
3. 调整粒子速度：将每个粒子的速度 𝑣𝑖乘以缩放因子 𝜆，得到新的速度 𝑣𝑖′
4. 更新系统动能和温度：通过新的速度计算新的动能和温度，检查是否接近目标温度。

Nose-Hoover 控温算法

Nose-Hoover 控温算法基于Nose-Hoover链（Nose-Hoover chain），其基本思想是通过引入一个与系统耦合的虚拟热浴来控制系统的动能（温度）。该虚拟热浴通过一个附加的自由度（即热浴变量）来模拟系统与外界环境的能量交换。具体的公式见书P116-118。 优点：

1. 精确控制温度：
   • Nose-Hoover 控温算法能够将系统的温度精确地保持在设定值上，而不会引入明显的温度波动。
2. 保留动力学性质：
   • 该算法不会破坏系统的动力学性质，可以准确地模拟出系统的时间相关行为和输运性质（如扩散系数）。
3. 能量守恒：
   • 系统在 Nose-Hoover 控温算法下仍然能保持总能量的守恒，只是在系统能量和热浴能量之间进行交换。

缺点：

1. 参数敏感性：
   • 算法对参数（如耦合常数）的选择较为敏感，参数选择不当可能导致控温效果不佳或数值不稳定。
2. 数值稳定性问题：
   • 对于某些系统，Nose-Hoover 控温算法可能导致数值不稳定，尤其是当系统较大或者系统中存在剧烈的动态变化时。
3. 初始态依赖性：
   • 算法对初始状态较为敏感，初始温度和虚拟热库参数的设置需要仔细调节。
4. 复杂性较高：
   • 与简单的控温方法（如Berendsen Thermostat或Langevin Thermostat）相比，Nose-Hoover 控温算法实现和调试相对复杂。

算法步骤：

1. 实施控温相关的时间演化算符 e*iLT∆t/2。
2. 实施与速度-Verlet 积分算法对应的时间演化算符。
3. 再次实施控温相关的时间演化算符 eiLT∆t/2。

为了克服NH控温算法的不足，Martyna、 Klein 和 Tuckerman 将NH算法中的一个控温变量扩展为一组控温变量，称为一条控温链 (chain)。所以，这种控温算法常称为NHC控温算法。在该算法中，第0个控温变量相当于NH控温算法中的那个控温变量，而第m(m ≥ 1)个控温变量的作用是控制第m − 1个控温变量的温度。详细公式见书P123-125

与传统的 Nose-Hoover 恒温器相比，NHC 方法具有以下优点：

更高的效率: NHC 方法可以更有效地控制系统的温度，特别是在模拟大型系统时。

更好的稳定性: NHC 方法比传统的 Nose-Hoover 恒温器更稳定，不易产生漂移。

更宽的温度范围: NHC 方法可以用于模拟更宽的温度范围内的系统。

Langevin Thermostat控温算法

郎之万控温算法（Langevin Thermostat）是一种用于分子动力学模拟中的控温方法，通过引入随机力和摩擦力来模拟与热浴（热库）的相互作用，从而使系统在一定的温度下进行演化。运动方程如下： $$\frac{d}{dt}{r}_i=\frac{p_i}{m_i}$$ $$d{p}_i=F_i{d}_t-\gamma p_i{d}_t+\sqrt{2\gamma m_i{k}_B{T}_o{d}_t}{R}_i$$

算法步骤：

1. 初始化粒子的位置和速度初始化粒子的位置和速度。设置模拟参数，如摩擦系数 𝛾、时间步长 Δ𝑡、目标温度 𝑇等。
2. 计算每个粒子所受的保守力F，计算摩擦力 F_friction和随机力R(t)
3. 使用数值积分方法（如Verlet算法或其改进版本）更新粒子的速度和位置。
4. 重复模拟循环2-3.

优点:

1. 控温效果好：郎之万控温算法能有效地控制系统的温度，使其在期望的温度附近波动，尤其适用于需要维持稳定温度的长时间模拟。
2. 简单易实现：该算法相对简单，容易在已有的分子动力学程序中实现。
3. 适用于小系统：在较小的系统中，郎之万控温算法的控温效果较好，因为摩擦和随机力可以有效地与系统中较少的

缺点:

1. 引入额外的力：摩擦力和随机力的引入可能会影响系统的动力学行为，导致结果偏离真实的物理过程。
2. 参数选择敏感：摩擦系数 ( \gamma ) 和随机力的强度需要合理选择，否则可能导致温度波动过大或过小，影响模拟的准确性。
3. 影响速度自相关函数：由于引入了非保守力，系统的速度自相关函数可能受到影响，从而影响运输性质（如扩散系数）的计算结果。

Andersen Thermostat 控温算法

Andersen Thermostat通过随机选择粒子并重置其速度来模拟与热浴（热库）的相互作用。其基本思想是通过随机碰撞来控制系统的温度，使其接近设定的目标温度。

计算思路：

1. 初始化：

   • 初始化粒子的速度和位置。
   • 设置模拟参数，如时间步长、目标温度、碰撞频率等。

2. 常规分子动力学步骤：

   • 按照分子动力学的常规方法，计算每个粒子所受的保守力。
   • 使用数值积分方法（如Verlet算法）更新粒子的速度和位置。

3. 随机碰撞步骤：

   • 对于每个时间步长，以概率来决定每个粒子是否与热浴发生碰撞。
   • 如果发生碰撞，重置该粒子的速度，新的速度从Maxwell-Boltzmann分布中随机采样。

优点

1. 控温效果好：能有效地将系统温度控制在目标温度附近。
2. 简单易实现：随机重置粒子的速度非常简单，不需要复杂的计算。
3. 对不同系统的适用性强：适用于各种类型的分子动力学模拟。

缺点

1. 破坏动量守恒：随机重置粒子的速度会破坏系统的动量守恒。
2. 对动力学行为的影响：频繁的速度重置会影响系统的动力学行为，可能不适用于某些需要精确动力学行为的模拟。

三、进行不同算法之间的对比

控温算法的性质对比：

书中考虑具有8000个原子的硅晶体超胞，采用Tersoff势函数描述原子间相互作用。采用GPUMD程序，测试几个控温算法获得平衡的时间。具体细节见书P135

四、控温知识扩展

在分子动力学模拟中，控温是确保系统热力学性质合理的重要部分。除了Berendsen控温算法、Nose-Hoover控温算法、Langevin Thermostat控温算法和Andersen Thermostat控温算法，还有一些其他控温方法和相关知识点对研究人员学习和应用分子动力学模拟非常重要。以下是一些扩展的控温知识：

1. Velocity Rescaling Thermostat（速度重标定控温算法）

速度重标定控温算法通过按比例缩放系统所有粒子的速度来调整温度，使系统的动能达到目标温度。这是一种简单但有效的控温方法，特别适用于初始阶段的温度调整。

公式 $$v_i\to v_i\sqrt{\frac{T_0}{T}}$$

其中 ( T ) 为当前温度， ( T_0 ) 为目标温度。

优点

• 简单易实现
• 适用于初始阶段快速调整温度

缺点

• 不保留系统的动力学性质
• 温度波动较大

2. Parrinello-Rahman Thermostat（Parrinello-Rahman 控温算法）

Parrinello-Rahman Thermostat 通过引入一个与系统耦合的虚拟热库来控制温度，类似于Nose-Hoover Thermostat，但该方法可以与Parrinello-Rahman Barostat一起使用，进行等温等压（NPT）系综的模拟。

优点

• 能同时控制温度和压力
• 适用于需要模拟真实物理条件的系统

缺点

• 实现复杂
• 参数选择敏感

3. Martyna-Tuckerman-Klein (MTK) Thermostat

MTK Thermostat 是对Nose-Hoover Thermostat的一种改进，通过引入多个热库变量，提高了数值稳定性。该方法特别适用于长时间的动力学模拟。

优点

• 控温精确
• 数值稳定性好
• 保留系统动力学性质

缺点

• 实现复杂
• 参数选择敏感

4. Stochastic Velocity Rescaling (SVR) Thermostat（随机速度重标定控温算法）

SVR Thermostat 结合了速度重标定和随机噪声，通过在速度重标定过程中引入随机扰动来控制温度。这种方法能更好地模拟真实热力学性质。

公式 $$v_i\to v_i\sqrt{\frac{T_0}{T}}+\sqrt{(1-\frac{T_0}{T})kT/m}{R}_i$$

优点

• 控温精确
• 保留系统动力学性质
• 适用于长时间模拟

缺点

• 算法实现相对复杂

扩展知识点

1. 控温算法的选择 不同的控温算法适用于不同的模拟需求和系统特性。研究人员应根据具体的模拟目标、系统性质和计算资源选择合适的控温算法。例如：

• 对初始温度调节要求较高的模拟，可以选择Berendsen或速度重标定控温算法。
• 对动力学性质要求高的长时间模拟，可以选择Nose-Hoover或Langevin Thermostat。
• 对需要模拟真实物理条件的NPT系综模拟，可以选择Parrinello-Rahman Thermostat。

2. 温度波动的控制 在实际模拟中，系统的温度会有一定的波动。合理的温度波动范围取决于系统大小和具体模拟的性质。过大的温度波动可能表明控温算法不稳定或参数选择不当。

3. 热库耦合常数的选择 热库耦合常数（如Nose-Hoover和Langevin Thermostat中的参数）决定了系统与热库的交换速率。合适的耦合常数能够平衡控温效果和系统动力学性质。

4. 长时间模拟的数值稳定性 在长时间的模拟中，数值稳定性是一个重要问题。选择合适的时间步长和控温参数能够提高模拟的数值稳定性，避免温度漂移和系统发散。