3. 分子动力学：从 NVE 到 NVT 系综
3.1. 引言
分子动力学（MD）是研究原子尺度动力学行为的核心方法。在计算物理中，我们需要在特定的热力学条件下研究材料性质——这就引出了统计系综的概念。
本章将从最基础的数值积分问题出发，逐步构建对MD系综的理解：
第一部分：探讨为什么MD需要特殊的数值积分器——通过对比实验展示保辛性的重要性
第二部分：深入理解微正则系综（NVE）的物理本质及其局限性
第三部分：引入恒温器概念，从简单到复杂逐步介绍各种NVT实现方法
第四部分：详解Nosé-Hoover链方法，理解扩展系统方法的精妙设计
3.1.1. 学习目标
完成本章学习后，你将能够：
理解保辛积分器的物理意义：为什么Velocity-Verlet比Runge-Kutta更适合MD
掌握算符分离思想：如何通过Trotter分解构建复杂的积分方案
认识温度控制的必要性：NVE系综的温度涨落问题及其对性质计算的影响
比较不同恒温器方法：理解各种NVT方法的物理基础、优缺点和适用场景
应用到实际问题：选择合适的系综和参数进行材料性质计算
3.1.2. 相关API文档
本章涉及的核心类和方法：
3.2. 环境准备
3.2.1. 系统设置说明
为了平衡教学清晰性和计算效率，我们采用分级系统策略：
小系统（2×2×2）：32原子，用于基础概念演示和算法细节展示
中等系统（3×3×3）：108原子，用于恒温器对比实验
较大系统（4×4×4）：256原子，用于统计性质收敛性演示
3.2.2. 势函数选择
本章使用EAM-Al1势函数，这是一个经过充分验证的铝材料势函数（Mendelev et al., 2008）。该势函数能够准确描述铝的弹性、热力学和缺陷性质。
# 导入必要的库和模块
import numpy as np
import plotly.graph_objects as go
# ThermoElasticSim核心模块
from thermoelasticsim.core.crystalline_structures import CrystallineStructureBuilder
from thermoelasticsim.core.structure import Cell
# MD传播子和积分方案
from thermoelasticsim.md.propagators import (
ForcePropagator,
PositionPropagator,
VelocityPropagator,
)
from thermoelasticsim.md.schemes import NVEScheme
# 势函数
from thermoelasticsim.potentials.eam import EAMAl1Potential
from thermoelasticsim.utils.utils import KB_IN_EV
# 初始化EAM势函数（需要C++扩展支持）
potential = EAMAl1Potential()
# 构建不同尺寸的FCC铝晶胞
builder = CrystallineStructureBuilder()
a0 = 4.045 # 铝的晶格常数，单位：Å
# 创建测试系统
cell32 = builder.create_fcc("Al", a0, (2, 2, 2)) # 32原子系统
print(f"小系统：{len(cell32.atoms)}个原子，用于基础演示")
# 定义基础工具函数
def kinetic_energy(cell: Cell) -> float:
"""计算系统动能"""
return float(cell.calculate_kinetic_energy())
def potential_energy(cell: Cell) -> float:
"""计算系统势能"""
return float(potential.calculate_energy(cell))
def temperature(cell: Cell) -> float:
"""计算瞬时温度"""
return float(cell.calculate_temperature())
def assign_maxwell_velocities(cell: Cell, temperature_K: float, seed: int = 42):
"""按Maxwell-Boltzmann分布初始化速度
Parameters
----------
cell : Cell
目标晶胞
temperature_K : float
目标温度（K）
seed : int
随机数种子，保证可重复性
"""
rng = np.random.default_rng(seed)
for atom in cell.atoms:
# Maxwell-Boltzmann分布的标准差：σ = sqrt(kT/m)
sigma = np.sqrt(KB_IN_EV * temperature_K / atom.mass)
atom.velocity = rng.normal(0.0, sigma, 3)
# 移除质心运动，避免系统漂移
cell.remove_com_motion()
print("环境准备完成！")
INFO: 字体配置成功，优先使用: Arial Unicode MS
小系统：32个原子，用于基础演示
环境准备完成！
# Plotly 输出助手：同时提供 vnd.plotly(JSON) 与 text/html 两种 MIME，兼容 VS Code 与 Sphinx
from IPython.display import display
class _PlotlyBundle:
def __init__(self, fig):
self.fig = fig
def _repr_mimebundle_(self, include=None, exclude=None):
# JSON for VS Code Jupyter plotly renderer
data_json = self.fig.to_plotly_json()
# HTML for Sphinx myst-nb off-mode（不内联JS，conf.py 注入 CDN）
html = self.fig.to_html(include_plotlyjs=False, full_html=False)
return {
"application/vnd.plotly.v1+json": data_json,
"text/html": html,
}
def pshow(fig):
display(_PlotlyBundle(fig))
3.3. 第一部分：数值积分与保辛性
3.4. 牛顿方程的数值挑战
分子动力学的核心是求解牛顿运动方程组：
\[m_i \frac{d^2\mathbf{r}_i}{dt^2} = \mathbf{F}_i = -\nabla_i U(\{\mathbf{r}\})\]
其中\(U\)是多体势能函数。这是一个高维（3N维）的二阶常微分方程组，通常没有解析解，必须采用数值方法。
3.4.1. 为什么普通数值方法不够好？
你可能会想：既然Runge-Kutta方法精度很高（四阶），为什么MD不用它？答案涉及一个深刻的物理原理：相空间体积守恒（刘维尔定理）。
在哈密顿力学框架下，系统在相空间中的演化必须保持相体积不变。这个性质被称为辛结构（symplectic structure）。保持这个结构的积分器称为保辛积分器（symplectic integrator）。
3.4.2. 保辛性的物理意义
能量守恒：虽然数值积分不能精确保持能量，但保辛积分器的能量误差有界，不会系统性增长
长期稳定性：保辛积分器在长时间模拟中保持物理性质的正确统计分布
时间可逆性：保辛积分器通常是时间可逆的，符合牛顿力学的基本对称性
3.5. 积分方法对比实验
让我们通过实验直观地看到不同积分方法的差异。我们将实现并比较：
欧拉方法：最简单的一阶方法，不保辛
Velocity-Verlet：MD标准方法，二阶精度，保辛
Runge-Kutta 4阶（可选）：高精度但不保辛
关键观察点：总能量\(E = T + V\)的长期行为。
def run_verlet(base: Cell, dt: float, steps: int):
"""使用Velocity-Verlet方法进行MD模拟
这是MD的标准积分方法，具有二阶精度和保辛性质。
算法步骤：
1. v(t+dt/2) = v(t) + F(t)/m * dt/2 （速度半步）
2. r(t+dt) = r(t) + v(t+dt/2) * dt （位置整步）
3. 计算F(t+dt) （更新力）
4. v(t+dt) = v(t+dt/2) + F(t+dt)/m * dt/2 （速度半步）
"""
cell = base.copy()
scheme = NVEScheme() # 使用库中实现的Velocity-Verlet
# 关键：必须先计算初始力
potential.calculate_forces(cell)
ts, Eks, Eps = [], [], []
t = 0.0
for s in range(steps):
scheme.step(cell, potential, dt)
if s % 1 == 0: # 每步记录
ts.append(t + dt)
Eks.append(kinetic_energy(cell))
Eps.append(potential_energy(cell))
t += dt
return np.array(ts), np.array(Eks), np.array(Eps)
def run_euler(base: Cell, dt: float, steps: int):
"""使用欧拉方法进行MD模拟（教学对比用）
最简单的积分方法，但不保辛，能量会系统性漂移。
算法步骤：
1. a(t) = F(t)/m （计算加速度）
2. v(t+dt) = v(t) + a(t) * dt （更新速度）
3. r(t+dt) = r(t) + v(t+dt) * dt （更新位置）
注意：这里使用的是"半隐式"欧拉（用新速度更新位置），
比完全显式欧拉稍微稳定一些。
"""
cell = base.copy()
potential.calculate_forces(cell)
ts, Eks, Eps = [], [], []
t = 0.0
for s in range(steps):
# 欧拉方法的核心：简单的前向差分
for atom in cell.atoms:
# 更新速度（使用当前时刻的力）
atom.velocity += (atom.force / atom.mass) * dt
# 更新位置（使用新速度）
atom.position += atom.velocity * dt
# 应用周期性边界条件
for atom in cell.atoms:
atom.position = cell.apply_periodic_boundary(atom.position)
# 重新计算力
potential.calculate_forces(cell)
if s % 1 == 0:
ts.append(t + dt)
Eks.append(kinetic_energy(cell))
Eps.append(potential_energy(cell))
t += dt
return np.array(ts), np.array(Eks), np.array(Eps)
# 执行对比实验
print("开始积分方法对比实验...")
# 准备相同的初始条件
base = cell32.copy()
assign_maxwell_velocities(base, 300.0, seed=123)
potential.calculate_forces(base) # 预计算初始力
# 模拟参数
dt = 0.2 # fs，较小的时间步长
steps = 5000 # 总时长 1000 fs
# 运行不同方法
print("运行欧拉方法...")
res_euler = run_euler(base, dt, steps)
print("运行Velocity-Verlet方法...")
res_verlet = run_verlet(base, dt, steps)
# 计算能量漂移
def calculate_drift(ts, Ek, Ep):
Et = Ek + Ep
drift = Et - Et[0]
drift_percent = (Et[-1] - Et[0]) / abs(Et[0]) * 100
return drift, drift_percent
drift_euler, pct_euler = calculate_drift(*res_euler)
drift_verlet, pct_verlet = calculate_drift(*res_verlet)
print("\n能量漂移分析（1000 fs后）：")
print(f"欧拉方法：{pct_euler:.3f}%")
print(f"Velocity-Verlet：{pct_verlet:.6f}%")
开始积分方法对比实验...
运行欧拉方法...
运行Velocity-Verlet方法...
能量漂移分析（1000 fs后）：
欧拉方法：0.003%
Velocity-Verlet：-0.000154%
# 可视化能量演化对比
fig = go.Figure()
# 添加总能量漂移曲线
fig.add_trace(
go.Scatter(
x=res_euler[0],
y=drift_euler,
mode="lines",
name="欧拉方法",
line=dict(color="#d62728", width=2),
)
)
fig.add_trace(
go.Scatter(
x=res_verlet[0],
y=drift_verlet,
mode="lines",
name="Velocity-Verlet",
line=dict(color="#2ca02c", width=2),
)
)
fig.update_layout(
title="保辛性对比：总能量漂移 ΔE(t) = E(t) - E(0)",
xaxis_title="时间 (fs)",
yaxis_title="能量漂移 (eV)",
template="plotly_white",
height=400,
annotations=[
dict(
text="欧拉方法：能量漂移（数值耗散）",
x=500,
y=drift_euler[2500],
showarrow=True,
arrowhead=2,
),
dict(
text="Velocity-Verlet：能量围绕零点振荡（保辛）",
x=500,
y=0,
showarrow=True,
arrowhead=2,
),
],
)
pshow(fig)
print("\n物理解释：")
print("- 欧拉方法的能量单调增长表明存在数值耗散，破坏了能量守恒")
print("- Velocity-Verlet的能量在零附近小幅振荡，长期平均守恒")
print("- 这就是为什么MD必须使用保辛积分器！")
物理解释：
- 欧拉方法的能量单调增长表明存在数值耗散，破坏了能量守恒
- Velocity-Verlet的能量在零附近小幅振荡，长期平均守恒
- 这就是为什么MD必须使用保辛积分器！
3.6. Velocity-Verlet算法详解
让我们深入理解Velocity-Verlet算法的精妙设计。这个算法可以从多个角度理解：
3.6.1. 1. 算法步骤分解
Velocity-Verlet算法将一个时间步分为四个子步骤：
速度半步更新：\(\mathbf{v}(t+\frac{dt}{2}) = \mathbf{v}(t) + \frac{\mathbf{F}(t)}{m}\frac{dt}{2}\)
位置整步更新：\(\mathbf{r}(t+dt) = \mathbf{r}(t) + \mathbf{v}(t+\frac{dt}{2})dt\)
力重新计算：\(\mathbf{F}(t+dt) = -\nabla U(\mathbf{r}(t+dt))\)
速度半步更新：\(\mathbf{v}(t+dt) = \mathbf{v}(t+\frac{dt}{2}) + \frac{\mathbf{F}(t+dt)}{m}\frac{dt}{2}\)
3.6.2. 2. 对称性与时间可逆性
注意算法的对称结构：速度更新被分成两个相同的半步，夹着位置更新。这种对称性保证了时间可逆性——如果你反转所有速度并运行相同时间，系统会精确回到初始状态（忽略数值舍入误差）。
# Velocity-Verlet单步演示：展示中间状态
print("Velocity-Verlet算法单步分解演示")
print("=" * 50)
# 准备测试系统
cell_demo = cell32.copy()
assign_maxwell_velocities(cell_demo, 300.0, seed=321)
# 关键：计算初始力
potential.calculate_forces(cell_demo)
print("初始状态：")
E0 = kinetic_energy(cell_demo) + potential_energy(cell_demo)
T0 = temperature(cell_demo)
print(f" 总能量 E = {E0:.6f} eV")
print(f" 温度 T = {T0:.1f} K")
# 创建传播子（Propagator）
v_prop = VelocityPropagator() # 速度更新
r_prop = PositionPropagator() # 位置更新
f_prop = ForcePropagator(potential) # 力计算
dt = 0.5 # fs
print(f"\n执行Velocity-Verlet步骤 (dt = {dt} fs)：")
# 步骤1：速度半步
print("1. 速度半步更新 v(t) → v(t+dt/2)")
v_prop.propagate(cell_demo, dt / 2)
print(f" 中间动能 = {kinetic_energy(cell_demo):.6f} eV")
# 步骤2：位置整步
print("2. 位置整步更新 r(t) → r(t+dt)")
r_prop.propagate(cell_demo, dt)
print(f" 新势能 = {potential_energy(cell_demo):.6f} eV")
# 步骤3：重新计算力
print("3. 计算新位置的力 F(t+dt)")
f_prop.propagate(cell_demo, dt) # dt参数在这里不使用
avg_force = np.mean([np.linalg.norm(atom.force) for atom in cell_demo.atoms])
print(f" 平均力大小 = {avg_force:.6f} eV/Å")
# 步骤4：速度半步
print("4. 速度半步更新 v(t+dt/2) → v(t+dt)")
v_prop.propagate(cell_demo, dt / 2)
# 最终状态
E1 = kinetic_energy(cell_demo) + potential_energy(cell_demo)
T1 = temperature(cell_demo)
print("\n最终状态：")
print(f" 总能量 E = {E1:.6f} eV")
print(f" 温度 T = {T1:.1f} K")
print(f"\n能量变化：ΔE = {E1 - E0:+.3e} eV ({abs(E1 - E0) / abs(E0) * 100:.4f}%)")
print("\n说明：能量的微小变化来自数值舍入误差，")
print("但这个误差有界且不会系统性增长。")
Velocity-Verlet算法单步分解演示
==================================================
初始状态：
总能量 E = -111.534543 eV
温度 T = 260.1 K
执行Velocity-Verlet步骤 (dt = 0.5 fs)：
1. 速度半步更新 v(t) → v(t+dt/2)
中间动能 = 1.043387 eV
2. 位置整步更新 r(t) → r(t+dt)
新势能 = -112.582587 eV
3. 计算新位置的力 F(t+dt)
平均力大小 = 0.136477 eV/Å
4. 速度半步更新 v(t+dt/2) → v(t+dt)
最终状态：
总能量 E = -111.536398 eV
温度 T = 261.1 K
能量变化：ΔE = -1.855e-03 eV (0.0017%)
说明：能量的微小变化来自数值舍入误差，
但这个误差有界且不会系统性增长。
3.7. 算符分离与Trotter分解
3.7.1. 理论基础
Velocity-Verlet算法实际上是算符分离（operator splitting）方法的一个特例。在哈密顿力学框架下，时间演化由刘维尔算符描述：
\[i\mathcal{L} = \sum_i \left( \frac{\mathbf{p}_i}{m_i} \cdot \frac{\partial}{\partial \mathbf{r}_i} + \mathbf{F}_i \cdot \frac{\partial}{\partial \mathbf{p}_i} \right)\]
这个算符可以分解为两部分：
\(i\mathcal{L}_T\)：动能部分（更新位置）
\(i\mathcal{L}_V\)：势能部分（更新动量/速度）
3.7.2. Trotter分解
时间演化算符的精确形式是：\(e^{i\mathcal{L}t}\)
由于\(\mathcal{L}_T\)和\(\mathcal{L}_V\)不对易，我们使用Trotter分解：
\[e^{i\mathcal{L}\Delta t} \approx e^{i\mathcal{L}_V \Delta t/2} e^{i\mathcal{L}_T \Delta t} e^{i\mathcal{L}_V \Delta t/2} + O(\Delta t^3)\]
这正是Velocity-Verlet算法的数学基础！
3.7.3. ThermoElasticSim的Propagator设计
基于算符分离思想，ThermoElasticSim采用了模块化的Propagator设计：
每个Propagator负责一个基本的时间演化操作
通过组合不同的Propagator，可以构建复杂的积分方案
这种设计使得从NVE扩展到NVT、NPT变得自然而优雅
具体实现见：
thermoelasticsim.md.propagators- 传播子模块thermoelasticsim.md.schemes- 积分方案模块step()- NVE时间步进方法
# 演示Propagator架构：手动组装 vs 使用NVEScheme
print("Propagator架构演示：比较手动组装和NVEScheme")
print("=" * 50)
# 准备两个相同的初始系统
cell_init = cell32.copy()
assign_maxwell_velocities(cell_init, 300.0, seed=777)
cell_manual = cell_init.copy() # 手动组装版本
cell_scheme = cell_init.copy() # 使用NVEScheme版本
# 都需要初始力
potential.calculate_forces(cell_manual)
potential.calculate_forces(cell_scheme)
# 手动组装Propagator
print("手动组装Propagator执行Velocity-Verlet：")
f_prop = ForcePropagator(potential)
v_prop = VelocityPropagator()
r_prop = PositionPropagator()
dt = 0.5
E0_manual = kinetic_energy(cell_manual) + potential_energy(cell_manual)
# 手动执行Velocity-Verlet
v_prop.propagate(cell_manual, dt / 2) # v半步
r_prop.propagate(cell_manual, dt) # r整步
f_prop.propagate(cell_manual, dt) # 更新力
v_prop.propagate(cell_manual, dt / 2) # v半步
E1_manual = kinetic_energy(cell_manual) + potential_energy(cell_manual)
# 使用NVEScheme
print("\n使用NVEScheme执行相同操作：")
scheme = NVEScheme()
E0_scheme = kinetic_energy(cell_scheme) + potential_energy(cell_scheme)
scheme.step(cell_scheme, potential, dt)
E1_scheme = kinetic_energy(cell_scheme) + potential_energy(cell_scheme)
# 比较结果
print("\n结果比较：")
print(f"手动组装：ΔE = {E1_manual - E0_manual:+.3e} eV")
print(f"NVEScheme：ΔE = {E1_scheme - E0_scheme:+.3e} eV")
print(f"差异：{abs(E1_manual - E1_scheme):.3e} eV")
print("\n结论：两种方法完全等价！")
print("NVEScheme内部就是按照相同顺序调用这些Propagator。")
print("\n这种模块化设计的优势：")
print("1. 代码清晰：每个Propagator职责单一")
print("2. 易于扩展：添加恒温器只需插入新的Propagator")
print("3. 保证正确性：对称性和时间可逆性由组合顺序保证")
Propagator架构演示：比较手动组装和NVEScheme
==================================================
手动组装Propagator执行Velocity-Verlet：
使用NVEScheme执行相同操作：
结果比较：
手动组装：ΔE = -3.926e-03 eV
NVEScheme：ΔE = -3.926e-03 eV
差异：0.000e+00 eV
结论：两种方法完全等价！
NVEScheme内部就是按照相同顺序调用这些Propagator。
这种模块化设计的优势：
1. 代码清晰：每个Propagator职责单一
2. 易于扩展：添加恒温器只需插入新的Propagator
3. 保证正确性：对称性和时间可逆性由组合顺序保证
3.8. 第二部分：NVE系综深入分析
3.9. 微正则系综的物理本质
NVE系综，也称为微正则系综，是最基础的统计系综：
N：粒子数固定
V：体积固定
E：总能量固定（理想情况）
在NVE系综中，系统是孤立的，与外界没有能量交换。这对应于理想的绝热条件。
3.9.1. 温度的统计意义
在NVE系综中，温度不是控制参数而是输出量。瞬时温度通过动能定义：
\[T = \frac{2K}{N_f k_B} = \frac{1}{N_f k_B} \sum_i m_i v_i^2\]
其中\(N_f = 3N - 3\)是自由度数（去除质心运动）。
3.9.2. NVE的局限性
温度涨落：瞬时温度围绕平均值涨落，涨落幅度\(\propto 1/\sqrt{N}\)
初始条件敏感：最终温度取决于初始能量分配
实验对应困难：大多数实验在恒温条件下进行
让我们通过长时间模拟观察这些特性。
# NVE长时间演化：观察能量守恒和温度涨落
print("NVE系综长时间演化分析")
print("=" * 50)
# 使用较大系统以获得更好的统计性质
cell256 = builder.create_fcc("Al", a0, (4, 4, 4)) # 256原子
print(f"系统规模：{len(cell256.atoms)}个原子")
# 初始化速度
for atom in cell256.atoms:
sigma = np.sqrt(KB_IN_EV * 300.0 / atom.mass)
atom.velocity = np.random.normal(0.0, sigma, 3)
cell256.remove_com_motion()
def run_nve_analysis(cell_init: Cell, dt: float, steps: int, sample_every: int = 10):
"""运行NVE模拟并收集统计数据"""
cell = cell_init.copy()
scheme = NVEScheme()
potential.calculate_forces(cell)
ts, Eks, Eps, Ts = [], [], [], []
t = 0.0
for s in range(steps):
scheme.step(cell, potential, dt)
if (s + 1) % sample_every == 0:
ts.append(t + dt)
Ek = kinetic_energy(cell)
Ep = potential_energy(cell)
T = temperature(cell)
Eks.append(Ek)
Eps.append(Ep)
Ts.append(T)
t += dt
return np.array(ts), np.array(Eks), np.array(Eps), np.array(Ts)
# 比较不同时间步长（控制总时长一致）
print("\n测试不同时间步长的稳定性：")
dt_values = [0.5, 1.0, 2.0] # fs
colors = ["#1f77b4", "#ff7f0e", "#2ca02c"]
data_series = {}
total_time = 1000.0 # 固定总时长 1000 fs
for dt in dt_values:
steps = int(total_time / dt) # 根据dt计算步数，保证总时长一致
print(f" 运行 dt = {dt} fs, 步数 = {steps}...")
ts, Ek, Ep, Tser = run_nve_analysis(cell256, dt=dt, steps=steps, sample_every=10)
data_series[dt] = (ts, Ek, Ep, Tser)
# 计算能量漂移
Et = Ek + Ep
drift_rate = (Et[-1] - Et[0]) / (ts[-1] * abs(Et[0])) * 1e6 # ppm/ps
print(f" 能量漂移率：{drift_rate:.3f} ppm/ps")
print(f" 温度：{np.mean(Tser):.1f} ± {np.std(Tser):.1f} K")
NVE系综长时间演化分析
==================================================
系统规模：256个原子
测试不同时间步长的稳定性：
运行 dt = 0.5 fs, 步数 = 2000...
能量漂移率：0.000 ppm/ps
温度：141.6 ± 24.8 K
运行 dt = 1.0 fs, 步数 = 1000...
能量漂移率：0.002 ppm/ps
温度：141.2 ± 23.8 K
运行 dt = 2.0 fs, 步数 = 500...
能量漂移率：0.000 ppm/ps
温度：140.4 ± 21.8 K
# 可视化：能量守恒对比
fig = go.Figure()
for i, (dt, color) in enumerate(zip(dt_values, colors, strict=False)):
ts, Ek, Ep, _ = data_series[dt]
Et = Ek + Ep
fig.add_trace(
go.Scatter(
x=ts,
y=(Et - Et[0]) * 1000, # 转换为meV
mode="lines",
name=f"dt = {dt} fs",
line=dict(color=color, width=2),
)
)
fig.update_layout(
title="NVE能量守恒 vs 时间步长",
xaxis_title="时间 (fs)",
yaxis_title="能量漂移 (meV)",
template="plotly_white",
height=400,
legend=dict(x=0.02, y=0.98),
)
pshow(fig)
print("\n观察：")
print("1. 较小的时间步长（0.5 fs）能量守恒最好")
print("2. 即使dt=2.0 fs，能量漂移仍然很小（<1 meV/ps）")
print("3. 这证明了Velocity-Verlet的稳健性")
观察：
1. 较小的时间步长（0.5 fs）能量守恒最好
2. 即使dt=2.0 fs，能量漂移仍然很小（<1 meV/ps）
3. 这证明了Velocity-Verlet的稳健性
# 温度涨落分析
# 使用dt=1.0的数据进行详细分析
ts, Ek, Ep, Tser = data_series[1.0]
# 计算统计量
T_mean = np.mean(Tser)
T_std = np.std(Tser)
T_rel_fluct = T_std / T_mean
print("温度涨落统计分析：")
print(f" 平均温度：{T_mean:.2f} K")
print(f" 标准差：{T_std:.2f} K")
print(f" 相对涨落：{T_rel_fluct * 100:.2f}%")
print(f" 理论预期（1/√N）：{1 / np.sqrt(len(cell256.atoms)) * 100:.2f}%")
# 绘制温度时间序列和分布
fig = go.Figure()
# 子图1：温度时间序列
fig.add_trace(
go.Scatter(
x=ts, y=Tser, mode="lines", name="瞬时温度", line=dict(color="#1f77b4", width=1)
)
)
# 添加平均值线
fig.add_hline(
y=T_mean,
line_dash="dash",
line_color="red",
annotation_text=f"平均值: {T_mean:.1f} K",
)
fig.update_layout(
title="NVE系综中的温度涨落",
xaxis_title="时间 (fs)",
yaxis_title="温度 (K)",
template="plotly_white",
height=400,
)
pshow(fig)
print("\n物理意义：")
print("- 温度涨落是NVE系综的固有特性")
print("- 涨落幅度与系统大小成反比（热力学极限下趋于零）")
print("- 这种涨落会影响性质计算的精度")
print("- 这就是为什么需要恒温器！")
温度涨落统计分析：
平均温度：141.23 K
标准差：23.84 K
相对涨落：16.88%
理论预期（1/√N）：6.25%
物理意义：
- 温度涨落是NVE系综的固有特性
- 涨落幅度与系统大小成反比（热力学极限下趋于零）
- 这种涨落会影响性质计算的精度
- 这就是为什么需要恒温器！
3.10. 第三部分：恒温器方法
3.11. 为什么需要恒温器？
实际应用中，我们通常需要在特定温度下研究材料性质。NVE系综的温度涨落和初始条件敏感性使其不适合这类研究。恒温器（thermostat）通过与虚拟热浴耦合，实现温度控制。
3.12. 恒温器的分类
根据实现原理，恒温器可分为三类：
3.12.1. 1. 速度重缩放方法
Berendsen恒温器：简单有效，但不产生正确的正则分布
适用于系统预平衡，不适合精确的热力学采样
3.12.2. 2. 随机方法
Andersen恒温器：随机碰撞，破坏动力学连续性
Langevin恒温器：摩擦+随机力，基于布朗运动理论
3.12.3. 3. 扩展系统方法
Nosé-Hoover恒温器：确定性方法，引入额外自由度
Nosé-Hoover链：改进的遍历性，消除周期性振荡
让我们逐一探讨这些方法。
# 准备通用测试函数
def prepare_cell(supercell=(3, 3, 3), T0=300.0, seed=2025):
"""准备测试晶胞
Parameters
----------
supercell : tuple
超胞大小
T0 : float
初始温度(K)
seed : int
随机数种子
"""
c = builder.create_fcc("Al", a0, supercell)
rng = np.random.default_rng(seed)
for atom in c.atoms:
sigma = np.sqrt(KB_IN_EV * T0 / atom.mass)
atom.velocity = rng.normal(0.0, sigma, 3)
c.remove_com_motion()
return c
def run_thermostat(scheme, cell: Cell, dt: float, steps: int, sample_every: int = 10):
"""运行恒温器模拟
Parameters
----------
scheme : IntegrationScheme
积分方案（NVT）
cell : Cell
初始晶胞
dt : float
时间步长(fs)
steps : int
总步数
sample_every : int
采样间隔
"""
c = cell.copy()
# 初始化力（重要！）
try:
potential.calculate_forces(c)
except Exception:
pass
ts, Tser, Ekser, Epser = [], [], [], []
t = 0.0
for s in range(steps):
scheme.step(c, potential, dt)
if (s + 1) % sample_every == 0:
ts.append(t + dt)
Tser.append(temperature(c))
Ekser.append(kinetic_energy(c))
Epser.append(potential_energy(c))
t += dt
return np.array(ts), np.array(Tser), np.array(Ekser), np.array(Epser)
print("恒温器测试环境准备完成！")
恒温器测试环境准备完成！
3.13. Berendsen恒温器
Berendsen恒温器是最简单的温度控制方法。它通过速度重缩放使系统温度指数衰减到目标值：
\[\lambda = \sqrt{1 + \frac{\Delta t}{\tau_T}\left(\frac{T_{\text{target}}}{T_{\text{current}}} - 1\right)}\]
其中\(\tau_T\)是温度弛豫时间。
3.13.1. 优点
实现简单，计算成本低
温度收敛快速平滑
适合系统初始平衡
3.13.2. 缺点
不产生正确的正则分布
抑制了能量涨落
不适合计算热力学性质
实现类：BerendsenNVTScheme
创建函数：create_berendsen_nvt_scheme()
# Berendsen恒温器演示
from thermoelasticsim.md.schemes import create_berendsen_nvt_scheme
print("Berendsen恒温器演示")
print("=" * 50)
# 准备系统：初始温度200K，目标温度300K
cell_ber = prepare_cell((3, 3, 3), T0=200.0, seed=10)
print(f"系统：{len(cell_ber.atoms)}个原子")
print("初始温度：200 K")
print("目标温度：300 K")
# 创建Berendsen恒温器
# tau参数：弛豫时间，通常取50-100倍时间步长
scheme_ber = create_berendsen_nvt_scheme(
target_temperature=300.0,
tau=100.0, # fs
)
# 运行模拟
dt = 0.5 # fs
steps = 4000
sample_every = 20
print("\n运行参数：")
print(f" 时间步长：{dt} fs")
print(f" 总时长：{dt * steps} fs")
print(" 弛豫时间τ：100 fs")
t_ber, T_ber, Ek_ber, Ep_ber = run_thermostat(
scheme_ber, cell_ber, dt, steps, sample_every
)
# 统计后半段温度
T_ber_tail = T_ber[len(T_ber) // 2 :]
print(f"\n平衡后温度：{np.mean(T_ber_tail):.1f} ± {np.std(T_ber_tail):.1f} K")
print(f"与目标偏差：{abs(np.mean(T_ber_tail) - 300.0):.1f} K")
# 可视化温度演化
import plotly.graph_objects as go
fig = go.Figure()
fig.add_trace(
go.Scatter(
x=t_ber,
y=T_ber,
mode="lines",
name="Berendsen",
line=dict(color="#1f77b4", width=2),
)
)
fig.add_hline(
y=300, line_dash="dash", line_color="red", annotation_text="目标温度 300 K"
)
fig.update_layout(
title="Berendsen恒温器温度演化",
xaxis_title="时间 (fs)",
yaxis_title="温度 (K)",
template="plotly_white",
height=400,
)
pshow(fig)
Berendsen恒温器演示
==================================================
系统：108个原子
初始温度：200 K
目标温度：300 K
Berendsen恒温器初始化: T=300.0K, τ_T=100.0fs, mode=equilibration
运行参数：
时间步长：0.5 fs
总时长：2000.0 fs
弛豫时间τ：100 fs
Berendsen统计 (#0步): T=183.3K, 误差=38.91%, 平均缩放=1.000
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#500步): T=209.2K, 误差=30.26%, 平均缩放=1.002
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#1000步): T=339.6K, 误差=13.20%, 平均缩放=1.001
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#1500步): T=279.2K, 误差=6.92%, 平均缩放=1.001
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#2000步): T=323.2K, 误差=7.74%, 平均缩放=1.001
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#2500步): T=297.8K, 误差=0.73%, 平均缩放=1.001
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#3000步): T=292.2K, 误差=2.61%, 平均缩放=1.000
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#3500步): T=311.9K, 误差=3.97%, 平均缩放=1.000
限制率: 0.00% (上限:0, 下限:0)
平衡后温度：300.2 ± 20.3 K
与目标偏差：0.2 K
3.14. Andersen恒温器
Andersen恒温器通过随机碰撞实现温度控制。在每个时间步，每个原子都有一定概率与虚拟热浴发生碰撞，碰撞后速度从Maxwell-Boltzmann分布重新采样。
3.14.1. 碰撞概率
小系统（<50原子）：\(\nu \approx 0.003-0.01\) fs\(^{-1}\)
大系统（>100原子）：\(\nu \approx 0.05-0.08\) fs\(^{-1}\)
3.14.2. 特点
严格产生正则分布
破坏动力学连续性
适合平衡态性质计算，不适合动力学性质
实现类：AndersenNVTScheme
创建函数：create_andersen_nvt_scheme()
# Andersen恒温器演示
from thermoelasticsim.md.schemes import create_andersen_nvt_scheme
print("Andersen恒温器演示")
print("=" * 50)
# 准备系统
cell_and = prepare_cell((3, 3, 3), T0=200.0, seed=11)
print(f"系统：{len(cell_and.atoms)}个原子")
# 创建Andersen恒温器
# collision_frequency：碰撞频率
scheme_and = create_andersen_nvt_scheme(
target_temperature=300.0,
collision_frequency=0.02, # fs^-1，中等系统的典型值
)
print("\n参数设置：")
print(" 碰撞频率：0.02 fs^-1")
print(f" 平均碰撞间隔：{1 / 0.02:.0f} fs")
print(f" 每步期望碰撞数：{0.02 * 0.5 * len(cell_and.atoms):.1f}")
# 运行模拟
t_and, T_and, Ek_and, Ep_and = run_thermostat(
scheme_and, cell_and, dt, steps, sample_every
)
# 分析温度控制
T_and_tail = T_and[len(T_and) // 2 :]
print(f"\n平衡后温度：{np.mean(T_and_tail):.1f} ± {np.std(T_and_tail):.1f} K")
print("\n说明：")
print("- Andersen方法的温度涨落比Berendsen大")
print("- 这是因为它产生正确的正则分布")
print("- 随机碰撞破坏了速度自相关函数")
# 可视化Andersen恒温器温度演化
fig = go.Figure()
fig.add_trace(
go.Scatter(
x=t_and,
y=T_and,
mode="lines",
name="Andersen",
line=dict(color="#ff7f0e", width=2),
)
)
fig.add_hline(
y=300, line_dash="dash", line_color="red", annotation_text="目标温度 300 K"
)
fig.update_layout(
title="Andersen恒温器温度演化（随机碰撞）",
xaxis_title="时间 (fs)",
yaxis_title="温度 (K)",
template="plotly_white",
height=400,
)
pshow(fig)
Andersen恒温器演示
==================================================
系统：108个原子
参数设置：
碰撞频率：0.02 fs^-1
平均碰撞间隔：50 fs
每步期望碰撞数：1.1
Andersen恒温器初始化: 108原子系统, 选择频率=0.02000 fs⁻¹
Andersen统计 (#1000步): T=286.5K, 误差=4.48%
每步碰撞数: 实际=1.0560, 期望=1.0800, 比值=0.978
Andersen统计 (#2000步): T=280.7K, 误差=6.43%
每步碰撞数: 实际=1.0610, 期望=1.0800, 比值=0.982
Andersen统计 (#3000步): T=285.5K, 误差=4.85%
每步碰撞数: 实际=1.0720, 期望=1.0800, 比值=0.993
Andersen统计 (#4000步): T=325.6K, 误差=8.53%
每步碰撞数: 实际=1.0793, 期望=1.0800, 比值=0.999
平衡后温度：287.3 ± 22.7 K
说明：
- Andersen方法的温度涨落比Berendsen大
- 这是因为它产生正确的正则分布
- 随机碰撞破坏了速度自相关函数
3.15. Langevin恒温器
Langevin恒温器基于布朗运动理论，通过摩擦力和随机力实现温度控制：
\[m\frac{d\mathbf{v}}{dt} = \mathbf{F} - \gamma m\mathbf{v} + \mathbf{R}(t)\]
其中：
\(\gamma\)：摩擦系数
\(\mathbf{R}(t)\)：随机力，满足涨落-耗散定理
3.15.1. 涨落-耗散定理
随机力必须满足：
\(\langle \mathbf{R}(t) \rangle = 0\)（零均值）
\(\langle \mathbf{R}(t)\mathbf{R}(t') \rangle = 2\gamma m k_B T \delta(t-t')\)（白噪声）
这保证了系统达到正确的热平衡。
3.15.2. 参数选择
弱耦合：\(\gamma \approx 0.01 ps^{-1}\)
中等耦合：\(\gamma \approx 0.1-1.0 ps^{-1}\)
强耦合：\(\gamma > 10 ps^{-1}\)（过阻尼极限）
实现类：LangevinNVTScheme
创建函数：create_langevin_nvt_scheme()
理论背景：分子动力学与统计系综
# Langevin恒温器演示
from thermoelasticsim.md.schemes import create_langevin_nvt_scheme
print("Langevin恒温器演示")
print("=" * 50)
# 准备系统
cell_lan = prepare_cell((3, 3, 3), T0=200.0, seed=12)
print(f"系统：{len(cell_lan.atoms)}个原子")
# 创建Langevin恒温器
# friction参数：摩擦系数，决定耦合强度
gamma = 2.0 # ps^-1，中等耦合
scheme_lan = create_langevin_nvt_scheme(target_temperature=300.0, friction=gamma)
print("\n参数设置：")
print(f" 摩擦系数γ：{gamma} ps^-1")
print(f" 弛豫时间τ：{1 / gamma:.2f} ps")
print(" 耦合强度：中等")
# 运行模拟
t_lan, T_lan, Ek_lan, Ep_lan = run_thermostat(
scheme_lan, cell_lan, dt, steps, sample_every
)
# 分析结果
T_lan_tail = T_lan[len(T_lan) // 2 :]
print(f"\n平衡后温度：{np.mean(T_lan_tail):.1f} ± {np.std(T_lan_tail):.1f} K")
# 计算总能量变化（不守恒是正常的）
Et_lan = Ek_lan + Ep_lan
dE_total = Et_lan[-1] - Et_lan[0]
print(f"\n总能量变化：{dE_total:.3f} eV")
print("\n说明：")
print("- Langevin方法中总能量不守恒（与热浴交换能量）")
print("- 摩擦力移除能量，随机力注入能量")
print("- 涨落-耗散定理保证达到正确的热平衡")
# 可视化Langevin恒温器温度演化
fig = go.Figure()
fig.add_trace(
go.Scatter(
x=t_lan,
y=T_lan,
mode="lines",
name="Langevin",
line=dict(color="#2ca02c", width=2),
)
)
fig.add_hline(
y=300, line_dash="dash", line_color="red", annotation_text="目标温度 300 K"
)
fig.update_layout(
title="Langevin恒温器温度演化（布朗动力学）",
xaxis_title="时间 (fs)",
yaxis_title="温度 (K)",
template="plotly_white",
height=400,
)
pshow(fig)
# 绘制能量变化
fig2 = go.Figure()
Et_lan = Ek_lan + Ep_lan
fig2.add_trace(
go.Scatter(
x=t_lan,
y=Et_lan - Et_lan[0],
mode="lines",
name="总能量变化",
line=dict(color="#9467bd", width=2),
)
)
fig2.update_layout(
title="Langevin恒温器：总能量变化（与热浴交换）",
xaxis_title="时间 (fs)",
yaxis_title="ΔE (eV)",
template="plotly_white",
height=400,
)
pshow(fig2)
Langevin恒温器演示
==================================================
系统：108个原子
Langevin恒温器初始化: T=300.0K, γ=2.000 ps⁻¹
阻尼时间: τ_damp = 0.500 ps
预期特性: 中等耦合
参数设置：
摩擦系数γ：2.0 ps^-1
弛豫时间τ：0.50 ps
耦合强度：中等
Langevin统计 (#1000步): T=259.4K, 误差=13.52%
能量平衡: 摩擦做功=-4.722e-03 eV, 随机做功=2.798e-04 eV
Langevin统计 (#2000步): T=265.7K, 误差=11.45%
能量平衡: 摩擦做功=-7.011e-03 eV, 随机做功=-5.732e-04 eV
Langevin统计 (#3000步): T=303.0K, 误差=1.00%
能量平衡: 摩擦做功=-7.768e-03 eV, 随机做功=3.631e-04 eV
Langevin统计 (#4000步): T=267.8K, 误差=10.73%
能量平衡: 摩擦做功=-7.952e-03 eV, 随机做功=-8.056e-04 eV
平衡后温度：284.2 ± 19.7 K
总能量变化：4.939 eV
说明：
- Langevin方法中总能量不守恒（与热浴交换能量）
- 摩擦力移除能量，随机力注入能量
- 涨落-耗散定理保证达到正确的热平衡
3.16. 第四部分：Nosé-Hoover链方法
3.17. 扩展系统方法的思想
Nosé-Hoover方法采用了完全不同的思路：不是通过外部干预（重缩放、碰撞、摩擦），而是将热浴作为系统的一部分，构建扩展的哈密顿系统。
3.17.1. 扩展哈密顿量
引入热浴变量\(\zeta\)及其共轭动量\(p_\zeta\)，扩展哈密顿量为：
\[H' = \sum_i \frac{\mathbf{p}_i^2}{2m_i} + U(\{\mathbf{r}\}) + \frac{p_\zeta^2}{2Q} + N_f k_B T \zeta\]
其中\(Q\)是热浴"质量"参数。
3.17.2. 运动方程
扩展系统的运动方程：
\(\dot{\mathbf{r}}_i = \mathbf{p}_i/m_i\)
\(\dot{\mathbf{p}}_i = \mathbf{F}_i - \zeta \mathbf{p}_i\)（摩擦项）
\(\dot{\zeta} = p_\zeta/Q\)
\(\dot{p}_\zeta = \sum_i \mathbf{p}_i^2/m_i - N_f k_B T\)（反馈项）
3.17.3. 链式扩展的必要性
单个Nosé-Hoover变量存在遍历性问题：系统可能陷入周期轨道。Nosé-Hoover链通过级联多个热浴变量解决这个问题。
# 单变量Nosé-Hoover演示：观察周期性振荡
from thermoelasticsim.md.schemes import NoseHooverNVTScheme
print("Nosé-Hoover恒温器演示")
print("=" * 50)
# 首先演示单变量情况（链长=1）
print("\n1. 单变量Nosé-Hoover（M=1）：")
cell_nh1 = prepare_cell((3, 3, 3), T0=200.0, seed=21)
# 创建单变量NH恒温器
scheme_nh1 = NoseHooverNVTScheme(
target_temperature=300.0,
tdamp=50.0, # fs，温度阻尼时间
tchain=1, # 链长=1（单变量）
tloop=1, # 内部循环次数
)
print(" 温度阻尼时间τ：50 fs")
print(" 链长M：1（单变量）")
# 运行较长时间以观察振荡
t_nh1, T_nh1, _, _ = run_thermostat(scheme_nh1, cell_nh1, dt, 6000, sample_every=20)
# 分析温度振荡
T_nh1_tail = T_nh1[len(T_nh1) // 2 :]
oscillation_amplitude = np.max(T_nh1_tail) - np.min(T_nh1_tail)
print("\n结果：")
print(f" 平均温度：{np.mean(T_nh1_tail):.1f} K")
print(f" 温度振荡幅度：{oscillation_amplitude:.1f} K")
print(f" 标准差：{np.std(T_nh1_tail):.1f} K")
if oscillation_amplitude > 50:
print(" ⚠️ 检测到明显的周期性振荡！")
# 可视化单变量NH的周期性振荡
fig = go.Figure()
fig.add_trace(
go.Scatter(
x=t_nh1,
y=T_nh1,
mode="lines",
name="NH (M=1)",
line=dict(color="#d62728", width=1.5),
)
)
fig.add_hline(
y=300, line_dash="dash", line_color="red", annotation_text="目标温度 300 K"
)
fig.update_layout(
title="单变量Nosé-Hoover：周期性温度振荡问题",
xaxis_title="时间 (fs)",
yaxis_title="温度 (K)",
template="plotly_white",
height=400,
annotations=[
dict(
text="明显的周期性振荡",
x=3000,
y=np.max(T_nh1[len(T_nh1) // 2 :]),
showarrow=True,
arrowhead=2,
)
],
)
pshow(fig)
Nosé-Hoover恒温器演示
==================================================
1. 单变量Nosé-Hoover（M=1）：
温度阻尼时间τ：50 fs
链长M：1（单变量）
NHC初始化: N_atoms=108, N_f=321, T=300.0K, τ=50.0fs
Q参数: Q[0]=2.075e+04
结果：
平均温度：310.5 K
温度振荡幅度：522.5 K
标准差：135.4 K
⚠️ 检测到明显的周期性振荡！
3.18. Nosé-Hoover链算法
3.18.1. 链式耦合
Nosé-Hoover链引入\(M\)个热浴变量\((\zeta_1, \zeta_2, ..., \zeta_M)\)，形成级联结构：
第一个变量\(\zeta_1\)与粒子系统耦合
每个\(\zeta_j\)与下一个\(\zeta_{j+1}\)耦合
最后一个\(\zeta_M\)自由演化
3.18.2. 质量参数选择
链变量的"质量"参数：
\(Q_1 = N_f k_B T \tau^2\)（第一个变量）
\(Q_j = k_B T \tau^2\)（\(j > 1\)）
其中\(\tau\)是特征时间尺度，通常取50-100个时间步。
3.18.3. 算符分离实现
NHC的时间演化通过多重时间尺度积分实现，使用Suzuki-Yoshida高阶分解提高精度。
实现类：NoseHooverNVTScheme
链传播子：NoseHooverChainPropagator
# Nosé-Hoover链（M=3）演示
print("\n2. Nosé-Hoover链（M=3）：")
cell_nhc = prepare_cell((3, 3, 3), T0=200.0, seed=22)
# 创建NHC恒温器
scheme_nhc = NoseHooverNVTScheme(
target_temperature=300.0,
tdamp=50.0, # fs
tchain=3, # 链长=3（推荐值）
tloop=1,
)
print(" 温度阻尼时间τ：50 fs")
print(" 链长M：3（消除振荡）")
# 运行模拟
t_nhc, T_nhc, _, _ = run_thermostat(scheme_nhc, cell_nhc, dt, 6000, sample_every=20)
# 分析结果
T_nhc_tail = T_nhc[len(T_nhc) // 2 :]
oscillation_amplitude_chain = np.max(T_nhc_tail) - np.min(T_nhc_tail)
print("\n结果：")
print(f" 平均温度：{np.mean(T_nhc_tail):.1f} K")
print(f" 温度振荡幅度：{oscillation_amplitude_chain:.1f} K")
print(f" 标准差：{np.std(T_nhc_tail):.1f} K")
print("\n改进效果：")
print(f" 单变量振荡幅度：{oscillation_amplitude:.1f} K")
print(f" 链式振荡幅度：{oscillation_amplitude_chain:.1f} K")
print(
f" 振荡减少：{(1 - oscillation_amplitude_chain / oscillation_amplitude) * 100:.1f}%"
)
print("\n结论：链式扩展显著改善了温度控制的稳定性！")
# 可视化NHC链的改进效果
fig = go.Figure()
# 添加单变量NH作为对比
fig.add_trace(
go.Scatter(
x=t_nh1,
y=T_nh1,
mode="lines",
name="NH (M=1)",
line=dict(color="#d62728", width=1, dash="dot"),
opacity=0.6,
)
)
# 添加NHC链
fig.add_trace(
go.Scatter(
x=t_nhc,
y=T_nhc,
mode="lines",
name="NHC (M=3)",
line=dict(color="#17becf", width=2),
)
)
fig.add_hline(
y=300, line_dash="dash", line_color="red", annotation_text="目标温度 300 K"
)
fig.update_layout(
title="Nosé-Hoover链：消除周期性振荡",
xaxis_title="时间 (fs)",
yaxis_title="温度 (K)",
template="plotly_white",
height=400,
legend=dict(x=0.7, y=0.95),
)
pshow(fig)
2. Nosé-Hoover链（M=3）：
温度阻尼时间τ：50 fs
链长M：3（消除振荡）
NHC初始化: N_atoms=108, N_f=321, T=300.0K, τ=50.0fs
Q参数: Q[0]=2.075e+04, Q[1]=6.463e+01
结果：
平均温度：298.5 K
温度振荡幅度：178.5 K
标准差：31.4 K
改进效果：
单变量振荡幅度：522.5 K
链式振荡幅度：178.5 K
振荡减少：65.8%
结论：链式扩展显著改善了温度控制的稳定性！
# NHC参数敏感性分析
print("\nNHC参数敏感性分析")
print("=" * 50)
# 测试不同的τ值
tau_values = [10.0, 100.0, 1000.0] # fs
colors_tau = ["#e377c2", "#7f7f7f", "#bcbd22"]
fig = go.Figure()
for tau, color in zip(tau_values, colors_tau, strict=False):
print(f"\n测试 τ = {tau} fs...")
cell_test = prepare_cell((3, 3, 3), T0=200.0, seed=30)
scheme_test = NoseHooverNVTScheme(
target_temperature=300.0, tdamp=tau, tchain=3, tloop=1
)
t_test, T_test, _, _ = run_thermostat(
scheme_test, cell_test, dt, 6000, sample_every=20
)
# 添加到图表
fig.add_trace(
go.Scatter(
x=t_test,
y=T_test,
mode="lines",
name=f"τ = {tau} fs",
line=dict(color=color, width=1.5),
)
)
# 统计
T_tail = T_test[len(T_test) // 2 :]
print(f" 平均温度：{np.mean(T_tail):.1f} ± {np.std(T_tail):.1f} K")
# 添加目标温度线
fig.add_hline(
y=300, line_dash="dash", line_color="red", annotation_text="目标温度 300 K"
)
fig.update_layout(
title="NHC参数敏感性：不同τ值的影响",
xaxis_title="时间 (fs)",
yaxis_title="温度 (K)",
template="plotly_white",
height=400,
)
pshow(fig)
print("\n参数选择指南：")
print("- τ太小（10 fs）：耦合太强，可能引起振荡")
print("- τ适中（100 fs）：平衡控制和响应速度")
print("- τ太大（1000 fs）：耦合太弱，收敛缓慢")
print("- 推荐：τ = 50-100 × dt")
NHC参数敏感性分析
==================================================
测试 τ = 10.0 fs...
NHC初始化: N_atoms=108, N_f=321, T=300.0K, τ=10.0fs
Q参数: Q[0]=8.298e+02, Q[1]=2.585e+00
平均温度：298.6 ± 24.7 K
测试 τ = 100.0 fs...
NHC初始化: N_atoms=108, N_f=321, T=300.0K, τ=100.0fs
Q参数: Q[0]=8.298e+04, Q[1]=2.585e+02
平均温度：293.5 ± 31.8 K
测试 τ = 1000.0 fs...
NHC初始化: N_atoms=108, N_f=321, T=300.0K, τ=1000.0fs
Q参数: Q[0]=8.298e+06, Q[1]=2.585e+04
平均温度：170.1 ± 27.5 K
参数选择指南：
- τ太小（10 fs）：耦合太强，可能引起振荡
- τ适中（100 fs）：平衡控制和响应速度
- τ太大（1000 fs）：耦合太弱，收敛缓慢
- 推荐：τ = 50-100 × dt
3.19. 第五部分：综合对比与应用
3.20. 恒温器方法全面对比
让我们在相同条件下比较所有恒温器方法，理解它们的特点和适用场景。
# 恒温器综合对比
print("恒温器方法综合对比")
print("=" * 50)
# 准备相同的初始系统
cell_base = prepare_cell((3, 3, 3), T0=200.0, seed=30)
print(f"测试系统：{len(cell_base.atoms)}个原子")
print("初始温度：200 K")
print("目标温度：300 K\n")
# 定义要比较的恒温器
thermostats = {
"Berendsen (τ=100fs)": create_berendsen_nvt_scheme(300.0, tau=100.0),
"Andersen (ν=0.02)": create_andersen_nvt_scheme(300.0, collision_frequency=0.02),
"Langevin (γ=2.0)": create_langevin_nvt_scheme(300.0, friction=2.0),
"NHC (τ=50fs, M=3)": NoseHooverNVTScheme(300.0, tdamp=50.0, tchain=3, tloop=1),
}
# 运行所有恒温器
results = {}
dt = 0.5
steps = 6000
sample_every = 20
for name, scheme in thermostats.items():
print(f"运行 {name}...")
cell = cell_base.copy()
ts, Tser, _, _ = run_thermostat(scheme, cell, dt, steps, sample_every)
results[name] = (ts, Tser)
# 创建对比图
fig = go.Figure()
colors = ["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728"]
for (name, (ts, Tser)), color in zip(results.items(), colors, strict=False):
fig.add_trace(
go.Scatter(
x=ts, y=Tser, mode="lines", name=name, line=dict(color=color, width=1.5)
)
)
# 添加目标温度线
fig.add_hline(y=300, line_dash="dash", line_color="black", annotation_text="目标 300 K")
fig.update_layout(
title="恒温器方法对比：温度演化",
xaxis_title="时间 (fs)",
yaxis_title="温度 (K)",
template="plotly_white",
height=500,
)
pshow(fig)
# 统计分析
print("\n统计分析（后半段数据）：")
print("-" * 60)
print(f"{'方法':<25} {'平均温度':<15} {'标准差':<10} {'误差'}")
print("-" * 60)
for name, (ts, Tser) in results.items():
# 只分析后半段（平衡后）
n = len(Tser)
tail = Tser[n // 2 :]
mean_T = np.mean(tail)
std_T = np.std(tail, ddof=1)
error = abs(mean_T - 300.0) / 300.0 * 100
print(f"{name:<25} {mean_T:<15.1f} {std_T:<10.1f} {error:.2f}%")
恒温器方法综合对比
==================================================
测试系统：108个原子
初始温度：200 K
目标温度：300 K
Berendsen恒温器初始化: T=300.0K, τ_T=100.0fs, mode=equilibration
Langevin恒温器初始化: T=300.0K, γ=2.000 ps⁻¹
阻尼时间: τ_damp = 0.500 ps
预期特性: 中等耦合
运行 Berendsen (τ=100fs)...
Berendsen统计 (#0步): T=207.0K, 误差=31.01%, 平均缩放=1.000
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#500步): T=208.0K, 误差=30.66%, 平均缩放=1.002
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#1000步): T=302.6K, 误差=0.88%, 平均缩放=1.001
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#1500步): T=277.6K, 误差=7.48%, 平均缩放=1.001
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#2000步): T=322.6K, 误差=7.54%, 平均缩放=1.001
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#2500步): T=281.7K, 误差=6.09%, 平均缩放=1.001
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#3000步): T=323.0K, 误差=7.66%, 平均缩放=1.000
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#3500步): T=314.6K, 误差=4.86%, 平均缩放=1.000
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#4000步): T=268.9K, 误差=10.36%, 平均缩放=1.000
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#4500步): T=320.0K, 误差=6.68%, 平均缩放=1.000
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#5000步): T=307.1K, 误差=2.37%, 平均缩放=1.000
限制率: 0.00% (上限:0, 下限:0)
Berendsen统计 (#5500步): T=313.0K, 误差=4.33%, 平均缩放=1.000
限制率: 0.00% (上限:0, 下限:0)
运行 Andersen (ν=0.02)...
Andersen恒温器初始化: 108原子系统, 选择频率=0.02000 fs⁻¹
Andersen统计 (#1000步): T=331.9K, 误差=10.65%
每步碰撞数: 实际=1.0890, 期望=1.0800, 比值=1.008
Andersen统计 (#2000步): T=323.8K, 误差=7.95%
每步碰撞数: 实际=1.0980, 期望=1.0800, 比值=1.017
Andersen统计 (#3000步): T=288.2K, 误差=3.95%
每步碰撞数: 实际=1.1030, 期望=1.0800, 比值=1.021
Andersen统计 (#4000步): T=302.1K, 误差=0.71%
每步碰撞数: 实际=1.1020, 期望=1.0800, 比值=1.020
Andersen统计 (#5000步): T=290.4K, 误差=3.21%
每步碰撞数: 实际=1.0988, 期望=1.0800, 比值=1.017
Andersen统计 (#6000步): T=305.2K, 误差=1.74%
每步碰撞数: 实际=1.0970, 期望=1.0800, 比值=1.016
运行 Langevin (γ=2.0)...
Langevin统计 (#1000步): T=236.3K, 误差=21.23%
能量平衡: 摩擦做功=-4.944e-03 eV, 随机做功=8.516e-05 eV
Langevin统计 (#2000步): T=277.7K, 误差=7.44%
能量平衡: 摩擦做功=-7.014e-03 eV, 随机做功=-5.146e-04 eV
Langevin统计 (#3000步): T=284.4K, 误差=5.21%
能量平衡: 摩擦做功=-7.481e-03 eV, 随机做功=1.912e-05 eV
Langevin统计 (#4000步): T=299.2K, 误差=0.27%
能量平衡: 摩擦做功=-7.835e-03 eV, 随机做功=-3.357e-04 eV
Langevin统计 (#5000步): T=291.9K, 误差=2.69%
能量平衡: 摩擦做功=-8.580e-03 eV, 随机做功=5.476e-04 eV
Langevin统计 (#6000步): T=262.4K, 误差=12.53%
能量平衡: 摩擦做功=-7.829e-03 eV, 随机做功=-1.497e-03 eV
运行 NHC (τ=50fs, M=3)...
NHC初始化: N_atoms=108, N_f=321, T=300.0K, τ=50.0fs
Q参数: Q[0]=2.075e+04, Q[1]=6.463e+01
统计分析（后半段数据）：
------------------------------------------------------------
方法 平均温度 标准差 误差
------------------------------------------------------------
Berendsen (τ=100fs) 300.3 19.7 0.11%
Andersen (ν=0.02) 288.3 22.9 3.91%
Langevin (γ=2.0) 292.3 20.2 2.57%
NHC (τ=50fs, M=3) 302.1 34.5 0.70%
3.21. 恒温器选择指南
基于上述对比，我们可以总结出恒温器的选择策略：
3.21.1. 使用场景推荐
场景 |
推荐方法 |
理由 |
|---|---|---|
系统预平衡 |
Berendsen |
快速收敛，计算便宜 |
平衡态采样 |
Langevin/NHC |
产生正确的正则分布 |
动力学性质 |
NHC |
保持动力学连续性 |
扩散系数 |
NVE或弱耦合Langevin |
最小化人工阻尼 |
自由能计算 |
NHC |
严格的系综采样 |
3.21.2. 参数选择建议
Berendsen：\(\tau_T = 50-100 \times dt\)
Andersen：根据系统大小调整碰撞频率
Langevin：\(\gamma = 0.1-2.0 ps^{-1}\)（平衡采样和动力学的折中）
NHC：\(\tau = 50-100 \times dt\)，链长\(M = 3\)
# 实际应用示例：计算径向分布函数
print("实际应用：不同恒温器下的径向分布函数")
print("=" * 50)
def calculate_rdf(cell, rmax=10.0, nbins=100):
"""计算径向分布函数（简化版）"""
positions = np.array([atom.position for atom in cell.atoms])
N = len(positions)
L = cell.lattice_vectors[0, 0] # 使用正确的属性获取晶格常数
# 初始化
dr = rmax / nbins
r = np.linspace(dr / 2, rmax - dr / 2, nbins)
g = np.zeros(nbins)
# 计算所有原子对距离
for i in range(N - 1):
for j in range(i + 1, N):
# 最小镜像约定
diff = positions[j] - positions[i]
diff = diff - L * np.round(diff / L)
dist = np.linalg.norm(diff)
if dist < rmax:
bin_idx = int(dist / dr)
if bin_idx < nbins:
g[bin_idx] += 2 # 计数两次（i-j和j-i）
# 归一化
for i in range(nbins):
shell_volume = 4 * np.pi * r[i] ** 2 * dr
g[i] = g[i] / (N * shell_volume * (N / L**3))
return r, g
print("计算平衡态RDF...")
print("（这里只是演示概念，实际应用需要更长的平衡和采样）\n")
# 使用NHC获得平衡构型
cell_eq = prepare_cell((3, 3, 3), T0=300.0, seed=40)
scheme_eq = NoseHooverNVTScheme(300.0, tdamp=50.0, tchain=3, tloop=1)
# 平衡运行
print("平衡系统...")
for _ in range(1000):
scheme_eq.step(cell_eq, potential, 0.5)
# 计算RDF
r, g = calculate_rdf(cell_eq, rmax=8.0, nbins=50)
# 绘制RDF
fig = go.Figure()
fig.add_trace(go.Scatter(x=r, y=g, mode="lines", line=dict(color="#1f77b4", width=2)))
fig.update_layout(
title="铝的径向分布函数 g(r)",
xaxis_title="距离 r (Å)",
yaxis_title="g(r)",
template="plotly_white",
height=400,
)
pshow(fig)
# 找出第一个峰
first_peak_idx = np.argmax(g)
first_peak_r = r[first_peak_idx]
print(f"\n第一近邻距离：{first_peak_r:.2f} Å")
print(f"FCC理论值：{a0 / np.sqrt(2):.2f} Å")
print(f"偏差：{abs(first_peak_r - a0 / np.sqrt(2)):.3f} Å")
实际应用：不同恒温器下的径向分布函数
==================================================
计算平衡态RDF...
（这里只是演示概念，实际应用需要更长的平衡和采样）
平衡系统...
NHC初始化: N_atoms=108, N_f=321, T=300.0K, τ=50.0fs
Q参数: Q[0]=2.075e+04, Q[1]=6.463e+01
第一近邻距离：2.80 Å
FCC理论值：2.86 Å
偏差：0.060 Å
3.22. 本章总结
3.22.1. 关键概念回顾
保辛积分器
Velocity-Verlet算法保持相空间体积和时间可逆性
长期能量稳定性优于高阶但非保辛的方法
算符分离和Trotter分解提供了理论框架
NVE系综
微正则系综是MD的基础
温度涨落是固有特性
适合研究动力学过程但不适合热力学采样
恒温器方法
不同方法有不同的物理基础和适用场景
Berendsen快速但不严格
Langevin基于布朗运动理论
NHC提供严格的正则采样
算符分离架构
Propagator设计使算法模块化
便于从NVE扩展到NVT、NPT
保证了算法的正确性和可维护性
3.22.2. 实践要点
时间步长选择：0.5-1.0 fs对大多数系统合适
系统大小效应：更大的系统有更小的涨落
平衡判断：监测温度、能量等物理量的收敛
参数优化：根据具体问题调整恒温器参数
3.22.3. 下章预告
在下一章（03b），我们将学习：
NPT系综：如何同时控制温度和压力
MTK算法：最严格的等温等压方法
算符分离的高级应用：处理更复杂的耦合
实际应用：相变、状态方程、弹性常数计算
通过本章的学习，你已经掌握了MD模拟的核心概念和基本系综。这些知识是进行材料性质计算的基础。继续前进，探索更丰富的MD世界！
3.22.4. 延伸阅读
分子动力学引擎 (md) - MD模块完整API文档
分子动力学与统计系综 - MD和系综理论基础
分子动力学：NPT 系综与 MTK 算法 - 下一章：NPT系综与MTK算法