4. 分子动力学：NPT 系综与 MTK 算法
4.1. 引言
在实际计算物理研究中，许多实验都是在恒定温度和压力条件下进行的。例如，材料的相变、热膨胀系数测量、高压物理实验等。NPT（等温等压）系综允许系统的体积随压力波动而变化，更真实地模拟实验条件。
本章将深入探讨NPT系综的理论基础和数值实现，重点介绍MTK（Martyna-Tobias-Klein）算法——目前最严格的NPT积分方案之一。
4.2. 学习目标
通过本章学习，您将：
理解压力控制的物理意义
为什么需要压力控制
压力的微观定义与计算
体积涨落与状态方程
掌握压力耦合方法
Berendsen压力耦合（简单但不严格）
Parrinello-Rahman方法
MTK算法的优势
深入理解MTK算法
扩展哈密顿量形式
九步对称Trotter分解
矩阵指数传播技术
实践NPT模拟
参数选择策略
收敛性诊断
高压模拟技巧
4.3. 第一部分：NPT系综理论基础
4.3.1. 1.1 吉布斯系综与压力控制
NPT系综（也称为等温等压系综或吉布斯系综）的配分函数为：
\[
\Delta_{NPT} = \frac{1}{N!h^{3N}} \int_0^\infty dV e^{-\beta PV} \int d\mathbf{r}^N d\mathbf{p}^N e^{-\beta H(\mathbf{r}, \mathbf{p})}
\]
其中：
\(P\) 是外部压力
\(V\) 是系统体积（作为动力学变量）
\(\beta = 1/(k_B T)\) 是逆温度
4.3.2. 1.2 压力的微观定义
在MD中，瞬时压力通过维里定理计算：
\[
P = \frac{1}{3V} \left( 2K + \sum_{i<j} \mathbf{r}_{ij} \cdot \mathbf{f}_{ij} \right)
\]
或等价地，通过应力张量：
\[
P = \frac{1}{3} \text{Tr}(\boldsymbol{\sigma})
\]
其中应力张量 \(\boldsymbol{\sigma}\) 包含动能和维里贡献。
4.3.3. 1.3 体积涨落与压缩率
在NPT系综中，体积涨落与等温压缩率相关：
\[
\langle (\Delta V)^2 \rangle = k_B T \kappa_T \langle V \rangle
\]
其中 \(\kappa_T = -\frac{1}{V}\left(\frac{\partial V}{\partial P}\right)_T\) 是等温压缩率。
# 导入必要的库和模块
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# ThermoElasticSim核心模块
from thermoelasticsim.core.crystalline_structures import CrystallineStructureBuilder
from thermoelasticsim.core.structure import Cell
from thermoelasticsim.elastic.mechanics import StressCalculator
from thermoelasticsim.md.schemes import (
create_mtk_npt_scheme,
)
from thermoelasticsim.potentials.eam import EAMAl1Potential
from thermoelasticsim.utils.utils import EV_TO_GPA, KB_IN_EV
print("ThermoElasticSim NPT教程环境已准备就绪")
print(f"玻尔兹曼常数 kB = {KB_IN_EV:.6f} eV/K")
print(f"压力转换因子 = {EV_TO_GPA:.3f} GPa/(eV/ų)")
INFO: 字体配置成功，优先使用: Arial Unicode MS
ThermoElasticSim NPT教程环境已准备就绪
玻尔兹曼常数 kB = 0.000086 eV/K
压力转换因子 = 160.218 GPa/(eV/ų)
# 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))
def prepare_cell(supercell=(3, 3, 3), T0=300.0, seed=2025):
"""
准备铝FCC晶胞并初始化速度
参数:
supercell: 超胞尺寸 (nx, ny, nz)
T0: 初始温度 (K)
seed: 随机数种子
返回:
Cell对象，速度已按Maxwell-Boltzmann分布初始化
"""
# 创建FCC铝结构
builder = CrystallineStructureBuilder()
a0 = 4.045 # Al的晶格常数 (Å)
cell = builder.create_fcc("Al", a0, supercell)
# Maxwell-Boltzmann速度初始化
rng = np.random.default_rng(seed)
for atom in cell.atoms:
# 速度标准差: σ = sqrt(kB*T/m)
sigma = np.sqrt(KB_IN_EV * T0 / atom.mass)
atom.velocity = rng.normal(0.0, sigma, 3)
# 移除质心运动（保持系统动量为零）
cell.remove_com_motion()
return cell
def compute_pressure_gpa(cell: Cell, potential) -> float:
"""
计算系统的瞬时压力
使用应力张量的迹: P = Tr(σ)/3
参数:
cell: Cell对象
potential: 势函数对象
返回:
压力值 (GPa)
"""
# 使用StressCalculator计算完整应力张量
calc = StressCalculator()
stress_tensor = calc.compute_stress(cell, potential)
# 压力 = 应力张量的迹 / 3
pressure_ev = np.trace(stress_tensor) / 3.0
# 转换为GPa
return float(pressure_ev * EV_TO_GPA)
def analyze_trajectory(times, values, label="Property"):
"""
分析轨迹数据的统计性质
"""
mean = np.mean(values)
std = np.std(values)
print(f"{label}:")
print(f" 平均值: {mean:.4f}")
print(f" 标准差: {std:.4f}")
print(f" 最小值: {np.min(values):.4f}")
print(f" 最大值: {np.max(values):.4f}")
return mean, std
# 初始化势函数和晶胞
pot = EAMAl1Potential()
print(f"使用势函数: {pot.__class__.__name__}")
# 创建3×3×3超胞（108个原子）
cell = prepare_cell((3, 3, 3), T0=300.0, seed=7)
print("\n系统信息:")
print(f" 原子数: {len(cell.atoms)}")
print(f" 初始体积: {cell.volume:.2f} ų")
print(" 晶格类型: FCC")
# 计算初始力和压力
pot.calculate_forces(cell)
P0 = compute_pressure_gpa(cell, pot)
T0 = cell.calculate_temperature()
print("\n初始状态:")
print(f" 温度: {T0:.2f} K")
print(f" 压力: {P0:+.4f} GPa")
print(f" 势能: {pot.calculate_energy(cell):.4f} eV")
print(f" 动能: {cell.calculate_kinetic_energy():.4f} eV")
使用势函数: EAMAl1Potential
系统信息:
原子数: 108
初始体积: 1786.98 ų
晶格类型: FCC
初始状态:
温度: 252.89 K
压力: -0.1217 GPa
势能: -368.3659 eV
动能: 3.4976 eV
4.4. 第二部分：压力耦合方法
4.4.1. 2.1 Berendsen压力耦合（教学演示）
Berendsen压力耦合是最简单的压力控制方法，通过缩放晶胞和原子坐标来调节压力：
\[
\mathbf{h} \rightarrow \mu \mathbf{h}, \quad \mathbf{r}_i \rightarrow \mu \mathbf{r}_i
\]
其中缩放因子：
\[
\mu = \left[1 - \kappa \frac{\Delta t}{\tau_P}(P - P_0)\right]^{1/3}
\]
\(\kappa\) 是压缩率（可调参数）
\(\tau_P\) 是压力弛豫时间
\(P_0\) 是目标压力
注意：Berendsen方法不产生正确的压力涨落，仅用于教学演示。实际应用请使用MTK算法。
def berendsen_pressure_coupling(
cell: Cell, potential, P_target: float, dt: float, tau_p: float, kappa: float = 5e-4
):
"""
Berendsen各向同性压力耦合（教学版）
参数:
cell: Cell对象
potential: 势函数
P_target: 目标压力 (GPa)
dt: 时间步长 (fs)
tau_p: 压力弛豫时间 (fs)
kappa: 有效压缩率 (1/GPa)
返回:
(当前压力, 缩放因子)
"""
# 计算当前压力
P_current = compute_pressure_gpa(cell, potential)
# 计算体积缩放因子
# μ³ = 1 - κ(P-P₀)Δt/τ_P
mu_cubed = 1.0 - kappa * (P_current - P_target) * (dt / tau_p)
mu = mu_cubed ** (1.0 / 3.0)
# 缩放晶格矢量
cell.lattice_vectors *= mu
# 缩放原子坐标（保持分数坐标不变）
for atom in cell.atoms:
atom.position *= mu
# 更新体积和力
cell.volume = cell.calculate_volume()
potential.calculate_forces(cell)
return P_current, mu
# 演示Berendsen压力耦合
print("演示Berendsen压力耦合:")
print("=" * 50)
# 准备新的晶胞
cell_ber = prepare_cell((3, 3, 3), T0=300.0, seed=10)
pot.calculate_forces(cell_ber)
# 模拟参数
dt = 0.5 # 时间步长 (fs)
tau_p = 5.0 # 压力弛豫时间 (fs)
P_target = 1.0 # 目标压力 (GPa)
steps = 2000 # 总步数
# 数据记录
times = []
pressures = []
volumes = []
# 运行Berendsen压力耦合
t = 0.0
for step in range(steps):
# 压力耦合步
P_now, scale = berendsen_pressure_coupling(cell_ber, pot, P_target, dt, tau_p)
# 记录数据（每10步）
if (step + 1) % 10 == 0:
times.append(t + dt)
pressures.append(P_now)
volumes.append(cell_ber.volume)
t += dt
# 分析结果
print(f"\n目标压力: {P_target:.2f} GPa")
print(f"最终压力: {pressures[-1]:.4f} GPa")
print(f"压力平均值: {np.mean(pressures[-100:]):.4f} GPa")
print(f"\n初始体积: {volumes[0]:.2f} ų")
print(f"最终体积: {volumes[-1]:.2f} ų")
print(f"体积变化: {(volumes[-1] / volumes[0] - 1) * 100:.2f}%")
演示Berendsen压力耦合:
==================================================
目标压力: 1.00 GPa
最终压力: 0.9988 GPa
压力平均值: 0.9897 GPa
初始体积: 1787.98 ų
最终体积: 1816.42 ų
体积变化: 1.59%
# 创建Berendsen结果的可视化
fig = make_subplots(
rows=2, cols=1, subplot_titles=("压力演化", "体积演化"), vertical_spacing=0.12
)
# 压力演化
fig.add_trace(
go.Scatter(
x=times, y=pressures, mode="lines", name="P(t)", line=dict(color="blue")
),
row=1,
col=1,
)
fig.add_hline(
y=P_target,
line_dash="dash",
line_color="red",
annotation_text=f"目标: {P_target} GPa",
row=1,
col=1,
)
# 体积演化
fig.add_trace(
go.Scatter(x=times, y=volumes, mode="lines", name="V(t)", line=dict(color="green")),
row=2,
col=1,
)
fig.update_xaxes(title_text="时间 (fs)", row=2, col=1)
fig.update_yaxes(title_text="压力 (GPa)", row=1, col=1)
fig.update_yaxes(title_text="体积 (ų)", row=2, col=1)
fig.update_layout(title="Berendsen压力耦合演示", height=600, showlegend=True)
pshow(fig)
4.5. 第三部分：MTK算法详解
4.5.1. 3.1 MTK扩展哈密顿量
MTK算法引入扩展相空间，包含粒子、晶胞和热浴/压浴变量：
\[
H_{MTK} = H_{particles} + H_{cell} + H_{thermostat} + H_{barostat}
\]
具体形式：
\[
H_{MTK} = \sum_i \frac{\mathbf{p}_i^2}{2m_i} + U(\mathbf{r}) + \frac{W_g \mathbf{p}_g^2}{2} + PV + H_{chains}
\]
其中：
\(\mathbf{p}_g\) 是晶胞动量（9个分量）
\(W_g\) 是晶胞质量参数
\(H_{chains}\) 是Nosé-Hoover链的贡献
4.5.2. 3.2 九步对称Trotter分解
MTK算法使用对称的九步分解保证时间可逆性：
压浴链传播 (Δt/2): 更新压力链变量
热浴链传播 (Δt/2): 更新温度链变量
晶胞动量更新 (Δt/2): \(\mathbf{p}_g \leftarrow \mathbf{p}_g + \frac{\Delta t}{2} \mathbf{F}_g\)
粒子动量更新 (Δt/2): 考虑晶胞耦合
位置和晶胞传播 (Δt): 矩阵指数更新
粒子动量更新 (Δt/2): 第二半步
晶胞动量更新 (Δt/2): 第二半步
热浴链传播 (Δt/2): 第二半步
压浴链传播 (Δt/2): 第二半步
这种对称分解确保了辛结构和长时间稳定性。
4.5.3. 3.3 参数选择指南
关键参数及推荐值：
参数 |
符号 |
推荐范围 |
说明 |
|---|---|---|---|
时间步长 |
Δt |
0.1-0.5 fs |
取决于势函数刚度 |
温度弛豫时间 |
τ_T |
20-100 fs |
越大温度涨落越小 |
压力弛豫时间 |
τ_P |
500-2000 fs |
越大压力涨落越小 |
热浴链长度 |
M_T |
3-5 |
通常3就足够 |
压浴链长度 |
M_P |
3-5 |
通常3就足够 |
# 使用 {py:class}`~thermoelasticsim.md.schemes.MTKNPTScheme` 进行NPT模拟
from thermoelasticsim.md.propagators import ForcePropagator
print("MTK算法九步分解演示")
print("=" * 50)
# 准备系统
cell_mtk = prepare_cell((3, 3, 3), T0=300.0, seed=12)
pot.calculate_forces(cell_mtk)
# 创建MTK NPT方案
# 使用 {py:func}`~thermoelasticsim.md.schemes.create_mtk_npt_scheme` 工厂函数
scheme = create_mtk_npt_scheme(
target_temperature=300.0, # 目标温度 (K)
target_pressure=0.0, # 目标压力 (GPa)
tdamp=50.0, # 温度弛豫时间 (fs)
pdamp=800.0, # 压力弛豫时间 (fs)
tchain=3, # 热浴链长度
pchain=3, # 压浴链长度
)
print("\nMTK参数设置:")
print(" 目标温度: 300.0 K")
print(" 目标压力: 0.0 GPa")
print(" 温度弛豫时间: 50.0 fs")
print(" 压力弛豫时间: 800.0 fs")
print(" 链长度: tchain=3, pchain=3")
# 演示单步九算符分解
dt = 0.2 # 时间步长
# 初始化内部传播子（通常自动完成）
if not hasattr(scheme, "_thermostat") or scheme._thermostat is None:
scheme._initialize_propagators(cell_mtk)
# 确保力传播子可用
if not hasattr(scheme, "_force_prop") or scheme._force_prop is None:
scheme._force_prop = ForcePropagator(pot)
# 记录每个子步骤的状态
def record_state(label, cell, pot):
"""记录系统状态"""
state = scheme.get_current_state(cell, pot)
return {
"step": label,
"T": state["temperature"],
"P": state["pressure"],
"V": state["volume"],
"E": state.get("conserved_energy", 0.0),
}
# 执行九步分解并记录
history = []
h = dt / 2 # 半步长
# 初始状态
history.append(record_state("0: 初始", cell_mtk, pot))
# 步骤1: 压浴链(半步)
scheme._barostat._integrate_barostat_chain(cell_mtk, h)
history.append(record_state("1: 压浴链½", cell_mtk, pot))
# 步骤2: 热浴链(半步)
scheme._thermostat.propagate(cell_mtk, h)
history.append(record_state("2: 热浴链½", cell_mtk, pot))
# 步骤3: 晶胞动量(半步)
scheme._barostat._update_cell_momenta(cell_mtk, pot, h)
history.append(record_state("3: p_cell½", cell_mtk, pot))
# 步骤4: 粒子动量(半步)
scheme._barostat.integrate_momenta(cell_mtk, pot, h)
history.append(record_state("4: p_atoms½", cell_mtk, pot))
# 步骤5: 位置和晶胞(全步)
scheme._barostat.propagate_positions_and_cell(cell_mtk, dt)
history.append(record_state("5: r+cell", cell_mtk, pot))
# 步骤6-9: 反向半步（保持对称性）
scheme._barostat.integrate_momenta(cell_mtk, pot, h)
history.append(record_state("6: p_atoms½", cell_mtk, pot))
scheme._barostat._update_cell_momenta(cell_mtk, pot, h)
history.append(record_state("7: p_cell½", cell_mtk, pot))
scheme._thermostat.propagate(cell_mtk, h)
history.append(record_state("8: 热浴链½", cell_mtk, pot))
scheme._barostat._integrate_barostat_chain(cell_mtk, h)
history.append(record_state("9: 压浴链½", cell_mtk, pot))
# 打印状态变化
print("\n九步分解状态演化:")
print(f"{'步骤':<15} {'温度(K)':<10} {'压力(GPa)':<12} {'体积(ų)':<10}")
print("-" * 50)
for h in history:
print(f"{h['step']:<15} {h['T']:<10.2f} {h['P']:<12.4f} {h['V']:<10.2f}")
MTK算法九步分解演示
==================================================
MTK-NPT积分方案初始化:
目标温度: 300.0 K
目标压力: 0.000 GPa
温度阻尼: 50.0 fs
压力阻尼: 800.0 fs
恒温器链: 3, 恒压器链: 3
MTK参数设置:
目标温度: 300.0 K
目标压力: 0.0 GPa
温度弛豫时间: 50.0 fs
压力弛豫时间: 800.0 fs
链长度: tchain=3, pchain=3
MTK恒压器初始化:
目标温度: 300.0 K
目标压力: 0.000000 eV/ų
压力阻尼: 800.0 fs
恒压器链: 3 变量
MTK-NPT传播子初始化完成
NHC初始化: N_atoms=108, N_f=321, T=300.0K, τ=50.0fs
Q参数: Q[0]=2.075e+04, Q[1]=6.463e+01
MTK参数初始化完成:
晶格质量 W: 1.80e+06
链质量 R[0]: 1.49e+05
链质量 R[1:]: 1.65e+04
九步分解状态演化:
步骤 温度(K) 压力(GPa) 体积(ų)
--------------------------------------------------
0: 初始 287.71 0.0000 1786.98
1: 压浴链½ 287.71 0.1505 1786.98
2: 热浴链½ 287.71 0.1505 1786.98
3: p_cell½ 287.71 0.1505 1786.98
4: p_atoms½ 287.71 0.1505 1786.98
5: r+cell 287.71 0.1523 1786.98
6: p_atoms½ 287.68 0.1523 1786.98
7: p_cell½ 287.68 0.1523 1786.98
8: 热浴链½ 287.68 0.1523 1786.98
9: 压浴链½ 287.68 0.1523 1786.98
# 可视化九步分解过程
steps = [h["step"] for h in history]
# 创建三个子图
fig = make_subplots(
rows=3,
cols=1,
subplot_titles=("温度变化", "压力变化", "体积变化"),
vertical_spacing=0.1,
)
# 温度
fig.add_trace(
go.Scatter(
x=steps,
y=[h["T"] for h in history],
mode="lines+markers",
name="T (K)",
line=dict(color="red"),
marker=dict(size=8),
),
row=1,
col=1,
)
# 压力
fig.add_trace(
go.Scatter(
x=steps,
y=[h["P"] for h in history],
mode="lines+markers",
name="P (GPa)",
line=dict(color="blue"),
marker=dict(size=8),
),
row=2,
col=1,
)
# 体积
fig.add_trace(
go.Scatter(
x=steps,
y=[h["V"] for h in history],
mode="lines+markers",
name="V (ų)",
line=dict(color="green"),
marker=dict(size=8),
),
row=3,
col=1,
)
fig.update_xaxes(title_text="MTK子步骤", row=3, col=1)
fig.update_yaxes(title_text="T (K)", row=1, col=1)
fig.update_yaxes(title_text="P (GPa)", row=2, col=1)
fig.update_yaxes(title_text="V (ų)", row=3, col=1)
fig.update_layout(title="MTK九步算符分解演示", height=800, showlegend=False)
pshow(fig)
4.5.4. 3.4 矩阵指数传播技术
MTK算法的核心技术之一是位置和晶胞的联合更新，通过矩阵指数实现：
\[\begin{split}
\begin{pmatrix}
\mathbf{r}(t+\Delta t) \\
\mathbf{h}(t+\Delta t)
\end{pmatrix} = \exp\left(\mathbf{G} \Delta t\right)
\begin{pmatrix}
\mathbf{r}(t) \\
\mathbf{h}(t)
\end{pmatrix}
\end{split}\]
其中 \(\mathbf{G}\) 是传播矩阵，包含晶胞速度信息。
矩阵指数的计算使用Padé近似或特征值分解，确保数值稳定性。
# 矩阵指数传播的简单演示（2D教学示例）
print("矩阵指数传播演示（2D简化版）")
print("=" * 50)
# 创建一个简单的2×2传播矩阵（旋转生成元）
A = np.array([[0.0, -1.0], [1.0, 0.0]]) # 逆时针旋转
dt = 0.2
print("\n传播矩阵 A:")
print(A)
# 计算矩阵指数 exp(A*dt)
from scipy.linalg import expm
expA = expm(A * dt)
print(f"\n矩阵指数 exp(A*{dt}):")
print(np.round(expA, 4))
# 应用到初始状态
x0 = np.array([1.0, 0.0]) # 初始位置
x1 = expA @ x0 # 传播后位置
print(f"\n初始状态 x(0) = {x0}")
print(f"传播后 x({dt}) = {np.round(x1, 4)}")
# 验证：这应该对应于角度 θ = dt 的旋转
theta = dt
x1_expected = np.array([np.cos(theta), np.sin(theta)])
print(f"理论值 = {np.round(x1_expected, 4)}")
print(f"误差 = {np.linalg.norm(x1 - x1_expected):.2e}")
# 可视化轨迹
n_steps = 50
trajectory = [x0]
x = x0.copy()
for _ in range(n_steps):
x = expA @ x
trajectory.append(x.copy())
trajectory = np.array(trajectory)
fig = go.Figure()
# 轨迹
fig.add_trace(
go.Scatter(
x=trajectory[:, 0],
y=trajectory[:, 1],
mode="lines+markers",
name="轨迹",
line=dict(color="blue"),
marker=dict(size=5),
)
)
# 起点和终点
fig.add_trace(
go.Scatter(
x=[trajectory[0, 0]],
y=[trajectory[0, 1]],
mode="markers",
name="起点",
marker=dict(size=10, color="green"),
)
)
fig.add_trace(
go.Scatter(
x=[trajectory[-1, 0]],
y=[trajectory[-1, 1]],
mode="markers",
name="终点",
marker=dict(size=10, color="red"),
)
)
fig.update_layout(
title="矩阵指数传播轨迹（2D旋转）",
xaxis_title="x",
yaxis_title="y",
xaxis=dict(scaleanchor="y", scaleratio=1),
width=600,
height=600,
)
pshow(fig)
矩阵指数传播演示（2D简化版）
==================================================
传播矩阵 A:
[[ 0. -1.]
[ 1. 0.]]
矩阵指数 exp(A*0.2):
[[ 0.9801 -0.1987]
[ 0.1987 0.9801]]
初始状态 x(0) = [1. 0.]
传播后 x(0.2) = [0.9801 0.1987]
理论值 = [0.9801 0.1987]
误差 = 2.78e-17
4.6. 第四部分：NPT模拟实践
4.6.1. 4.1 参数敏感性分析
压力弛豫时间 \(\tau_P\) 是NPT模拟的关键参数：
\(\tau_P\) 太小：压力振荡剧烈，难以收敛
\(\tau_P\) 太大：响应缓慢，需要更长模拟时间
推荐值：500-1000 fs（金属体系）
# 测试不同pdamp值的影响
def run_mtk_with_pdamp(pdamp_val: float, steps=3000, sample_every=20):
"""
运行MTK NPT模拟，测试pdamp参数影响
参数:
pdamp_val: 压力弛豫时间 (fs)
steps: 总步数
sample_every: 采样间隔
"""
# 准备新系统
c = prepare_cell((3, 3, 3), T0=300.0, seed=33)
pot.calculate_forces(c)
# 创建MTK方案
scheme = create_mtk_npt_scheme(
target_temperature=300.0,
target_pressure=0.0,
tdamp=50.0,
pdamp=pdamp_val, # 变化的参数
tchain=3,
pchain=3,
)
dt = 0.2
times = []
pressures = []
t = 0.0
# 运行模拟
for step in range(steps):
scheme.step(c, pot, dt)
if (step + 1) % sample_every == 0:
state = scheme.get_current_state(c, pot)
times.append(t + dt)
pressures.append(state["pressure"])
t += dt
return np.array(times), np.array(pressures)
print("测试不同pdamp值的影响...")
print("=" * 50)
# 测试三个不同的pdamp值
pdamp_values = [250.0, 500.0, 1000.0]
results = {}
for pdamp in pdamp_values:
print(f"\n运行 pdamp = {pdamp} fs...")
times, pressures = run_mtk_with_pdamp(pdamp)
results[pdamp] = (times, pressures)
# 计算统计
p_mean = np.mean(pressures[-50:])
p_std = np.std(pressures[-50:])
print(f" 最后50点: 平均压力 = {p_mean:.4f} ± {p_std:.4f} GPa")
# 可视化比较
fig = go.Figure()
colors = ["#1f77b4", "#ff7f0e", "#2ca02c"]
for (pdamp, (times, pressures)), color in zip(results.items(), colors, strict=False):
fig.add_trace(
go.Scatter(
x=times,
y=pressures,
mode="lines",
name=f"τ_P = {pdamp} fs",
line=dict(color=color, width=2),
)
)
# 添加目标压力线
fig.add_hline(y=0.0, line_dash="dash", line_color="red", annotation_text="目标: 0 GPa")
fig.update_layout(
title="MTK压力演化：τ_P敏感性分析",
xaxis_title="时间 (fs)",
yaxis_title="压力 (GPa)",
height=500,
legend=dict(x=0.7, y=1),
)
pshow(fig)
测试不同pdamp值的影响...
==================================================
运行 pdamp = 250.0 fs...
MTK-NPT积分方案初始化:
目标温度: 300.0 K
目标压力: 0.000 GPa
温度阻尼: 50.0 fs
压力阻尼: 250.0 fs
恒温器链: 3, 恒压器链: 3
MTK恒压器初始化:
目标温度: 300.0 K
目标压力: 0.000000 eV/ų
压力阻尼: 250.0 fs
恒压器链: 3 变量
MTK-NPT传播子初始化完成
MTK参数初始化完成:
晶格质量 W: 1.76e+05
链质量 R[0]: 1.45e+04
链质量 R[1:]: 1.62e+03
NHC初始化: N_atoms=108, N_f=321, T=300.0K, τ=50.0fs
Q参数: Q[0]=2.075e+04, Q[1]=6.463e+01
最后50点: 平均压力 = 0.0050 ± 0.1430 GPa
运行 pdamp = 500.0 fs...
MTK-NPT积分方案初始化:
目标温度: 300.0 K
目标压力: 0.000 GPa
温度阻尼: 50.0 fs
压力阻尼: 500.0 fs
恒温器链: 3, 恒压器链: 3
MTK恒压器初始化:
目标温度: 300.0 K
目标压力: 0.000000 eV/ų
压力阻尼: 500.0 fs
恒压器链: 3 变量
MTK-NPT传播子初始化完成
MTK参数初始化完成:
晶格质量 W: 7.04e+05
链质量 R[0]: 5.82e+04
链质量 R[1:]: 6.46e+03
NHC初始化: N_atoms=108, N_f=321, T=300.0K, τ=50.0fs
Q参数: Q[0]=2.075e+04, Q[1]=6.463e+01
最后50点: 平均压力 = 0.0319 ± 0.6069 GPa
运行 pdamp = 1000.0 fs...
MTK-NPT积分方案初始化:
目标温度: 300.0 K
目标压力: 0.000 GPa
温度阻尼: 50.0 fs
压力阻尼: 1000.0 fs
恒温器链: 3, 恒压器链: 3
MTK恒压器初始化:
目标温度: 300.0 K
目标压力: 0.000000 eV/ų
压力阻尼: 1000.0 fs
恒压器链: 3 变量
MTK-NPT传播子初始化完成
MTK参数初始化完成:
晶格质量 W: 2.82e+06
链质量 R[0]: 2.33e+05
链质量 R[1:]: 2.59e+04
NHC初始化: N_atoms=108, N_f=321, T=300.0K, τ=50.0fs
Q参数: Q[0]=2.075e+04, Q[1]=6.463e+01
最后50点: 平均压力 = 0.0412 ± 2.0596 GPa
4.6.2. 4.2 完整NPT模拟示例
现在进行一个完整的NPT模拟，监控所有重要物理量：
温度和压力收敛
体积演化
扩展守恒量（用于验证积分器正确性）
# 完整NPT模拟：0 GPa, 300 K
print("运行完整NPT模拟 (T=300K, P=0GPa)")
print("=" * 50)
# 准备系统
cell_npt = prepare_cell((3, 3, 3), T0=300.0, seed=40)
pot.calculate_forces(cell_npt)
# 创建MTK方案（优化参数）
scheme_npt = create_mtk_npt_scheme(
target_temperature=300.0,
target_pressure=0.0,
tdamp=50.0, # 适中的温度弛豫
pdamp=800.0, # 较大的压力弛豫（更平滑）
tchain=3,
pchain=3,
)
# 模拟参数
dt = 0.2
n_steps = 5000
sample_interval = 20
# 数据收集
trajectory = {
"time": [],
"temperature": [],
"pressure": [],
"volume": [],
"conserved": [],
"potential": [],
"kinetic": [],
}
# 运行模拟
t = 0.0
print("\n模拟进度:")
for step in range(n_steps):
# MTK积分步
scheme_npt.step(cell_npt, pot, dt)
# 定期采样
if (step + 1) % sample_interval == 0:
state = scheme_npt.get_current_state(cell_npt, pot)
trajectory["time"].append(t + dt)
trajectory["temperature"].append(state["temperature"])
trajectory["pressure"].append(state["pressure"])
trajectory["volume"].append(state["volume"])
trajectory["conserved"].append(state.get("conserved_energy", 0))
trajectory["potential"].append(pot.calculate_energy(cell_npt))
trajectory["kinetic"].append(cell_npt.calculate_kinetic_energy())
# 进度报告
if (step + 1) % 1000 == 0:
print(f" 步骤 {step + 1}/{n_steps} 完成")
t += dt
# 分析结果
print("\n模拟结果分析:")
print("=" * 50)
# 使用后半部分数据计算平衡性质
equil_start = len(trajectory["time"]) // 2
T_mean = np.mean(trajectory["temperature"][equil_start:])
T_std = np.std(trajectory["temperature"][equil_start:])
P_mean = np.mean(trajectory["pressure"][equil_start:])
P_std = np.std(trajectory["pressure"][equil_start:])
V_mean = np.mean(trajectory["volume"][equil_start:])
V_std = np.std(trajectory["volume"][equil_start:])
print("\n平衡态性质（后半程平均）:")
print(f" 温度: {T_mean:.2f} ± {T_std:.2f} K")
print(f" 压力: {P_mean:.4f} ± {P_std:.4f} GPa")
print(f" 体积: {V_mean:.2f} ± {V_std:.2f} ų")
# 检查守恒量漂移
if trajectory["conserved"][0] != 0:
drift = (trajectory["conserved"][-1] - trajectory["conserved"][0]) / trajectory[
"conserved"
][0]
print(f"\n扩展守恒量相对漂移: {drift:.2e}")
运行完整NPT模拟 (T=300K, P=0GPa)
==================================================
MTK-NPT积分方案初始化:
目标温度: 300.0 K
目标压力: 0.000 GPa
温度阻尼: 50.0 fs
压力阻尼: 800.0 fs
恒温器链: 3, 恒压器链: 3
模拟进度:
MTK恒压器初始化:
目标温度: 300.0 K
目标压力: 0.000000 eV/ų
压力阻尼: 800.0 fs
恒压器链: 3 变量
MTK-NPT传播子初始化完成
MTK参数初始化完成:
晶格质量 W: 1.80e+06
链质量 R[0]: 1.49e+05
链质量 R[1:]: 1.65e+04
NHC初始化: N_atoms=108, N_f=321, T=300.0K, τ=50.0fs
Q参数: Q[0]=2.075e+04, Q[1]=6.463e+01
步骤 1000/5000 完成
步骤 2000/5000 完成
步骤 3000/5000 完成
步骤 4000/5000 完成
步骤 5000/5000 完成
模拟结果分析:
==================================================
平衡态性质（后半程平均）:
温度: 299.08 ± 43.64 K
压力: 0.0251 ± 0.6006 GPa
体积: 1825.42 ± 11.29 ų
扩展守恒量相对漂移: -2.11e-02
# 创建综合可视化
fig = make_subplots(
rows=3,
cols=2,
subplot_titles=(
"温度演化",
"压力演化",
"体积演化",
"势能演化",
"动能演化",
"守恒量演化",
),
vertical_spacing=0.15,
horizontal_spacing=0.12,
)
# 温度
fig.add_trace(
go.Scatter(
x=trajectory["time"],
y=trajectory["temperature"],
mode="lines",
name="T(t)",
line=dict(color="red", width=1),
),
row=1,
col=1,
)
fig.add_hline(y=300, line_dash="dash", line_color="darkred", row=1, col=1)
# 压力
fig.add_trace(
go.Scatter(
x=trajectory["time"],
y=trajectory["pressure"],
mode="lines",
name="P(t)",
line=dict(color="blue", width=1),
),
row=1,
col=2,
)
fig.add_hline(y=0, line_dash="dash", line_color="darkblue", row=1, col=2)
# 体积
fig.add_trace(
go.Scatter(
x=trajectory["time"],
y=trajectory["volume"],
mode="lines",
name="V(t)",
line=dict(color="green", width=1),
),
row=2,
col=1,
)
# 势能
fig.add_trace(
go.Scatter(
x=trajectory["time"],
y=trajectory["potential"],
mode="lines",
name="U(t)",
line=dict(color="purple", width=1),
),
row=2,
col=2,
)
# 动能
fig.add_trace(
go.Scatter(
x=trajectory["time"],
y=trajectory["kinetic"],
mode="lines",
name="K(t)",
line=dict(color="orange", width=1),
),
row=3,
col=1,
)
# 守恒量
fig.add_trace(
go.Scatter(
x=trajectory["time"],
y=trajectory["conserved"],
mode="lines",
name="H(t)",
line=dict(color="black", width=1),
),
row=3,
col=2,
)
# 更新坐标轴标签
fig.update_xaxes(title_text="时间 (fs)", row=2, col=1)
fig.update_xaxes(title_text="时间 (fs)", row=2, col=2)
fig.update_yaxes(title_text="T (K)", row=1, col=1)
fig.update_yaxes(title_text="P (GPa)", row=1, col=2)
fig.update_yaxes(title_text="V (ų)", row=2, col=1)
fig.update_yaxes(title_text="能量 (eV)", row=2, col=2)
fig.update_layout(title="NPT模拟结果 (300K, 0GPa)", height=700, showlegend=True)
pshow(fig)
4.6.3. 4.3 高压NPT模拟
测试MTK算法在高压条件下的表现。高压会导致：
体积显著减小
原子间距离缩短
势能增加
可能需要更小的时间步长
# 高压NPT模拟: 10 GPa
print("高压NPT模拟 (T=300K, P=10GPa)")
print("=" * 50)
# 准备较大的系统以获得更好的统计
cell_hp = prepare_cell((4, 4, 4), T0=200.0, seed=41) # 256原子
pot.calculate_forces(cell_hp)
print(f"\n系统规模: {len(cell_hp.atoms)} 原子")
print(f"初始体积: {cell_hp.volume:.2f} ų")
# 创建高压MTK方案
scheme_hp = create_mtk_npt_scheme(
target_temperature=300.0,
target_pressure=10.0, # 10 GPa高压
tdamp=50.0,
pdamp=100.0, # 较小的pdamp以快速响应
tchain=3,
pchain=3,
)
# 模拟参数
dt = 0.2
n_steps = 8000
sample_interval = 40
# 数据收集
hp_data = {
"time": [],
"pressure": [],
"volume": [],
"density": [], # 添加密度跟踪
}
# 计算初始密度
total_mass = sum(atom.mass for atom in cell_hp.atoms)
initial_density = total_mass / cell_hp.volume # amu/ų
# 运行高压模拟
t = 0.0
print("\n运行高压模拟...")
for step in range(n_steps):
scheme_hp.step(cell_hp, pot, dt)
if (step + 1) % sample_interval == 0:
state = scheme_hp.get_current_state(cell_hp, pot)
density = total_mass / state["volume"]
hp_data["time"].append(t + dt)
hp_data["pressure"].append(state["pressure"])
hp_data["volume"].append(state["volume"])
hp_data["density"].append(density)
if (step + 1) % 1000 == 0:
current_p = scheme_hp.get_current_state(cell_hp, pot)["pressure"]
print(f" 步骤 {step + 1}: P = {current_p:.2f} GPa")
t += dt
# 分析高压结果
equil_start = len(hp_data["time"]) // 2
final_volume = np.mean(hp_data["volume"][equil_start:])
final_pressure = np.mean(hp_data["pressure"][equil_start:])
final_density = np.mean(hp_data["density"][equil_start:])
print("\n高压平衡态结果:")
print(f" 最终压力: {final_pressure:.2f} GPa")
print(f" 体积变化: {(final_volume / hp_data['volume'][0] - 1) * 100:.2f}%")
print(f" 密度变化: {(final_density / initial_density - 1) * 100:.2f}%")
# 估算体积模量
delta_P = final_pressure - 0.0 # 相对于0 GPa
delta_V_V = (final_volume - hp_data["volume"][0]) / hp_data["volume"][0]
if delta_V_V != 0:
bulk_modulus = delta_P / delta_V_V
print(f" 估算体积模量: {bulk_modulus:.1f} GPa")
高压NPT模拟 (T=300K, P=10GPa)
==================================================
系统规模: 256 原子
初始体积: 4235.80 ų
MTK-NPT积分方案初始化:
目标温度: 300.0 K
目标压力: 10.000 GPa
温度阻尼: 50.0 fs
压力阻尼: 100.0 fs
恒温器链: 3, 恒压器链: 3
运行高压模拟...
MTK恒压器初始化:
目标温度: 300.0 K
目标压力: 0.062415 eV/ų
压力阻尼: 100.0 fs
恒压器链: 3 变量
MTK-NPT传播子初始化完成
MTK参数初始化完成:
晶格质量 W: 6.64e+04
链质量 R[0]: 2.33e+03
链质量 R[1:]: 2.59e+02
NHC初始化: N_atoms=256, N_f=765, T=300.0K, τ=50.0fs
Q参数: Q[0]=4.944e+04, Q[1]=6.463e+01
步骤 1000: P = 15.50 GPa
步骤 2000: P = 1.77 GPa
步骤 3000: P = 19.91 GPa
步骤 4000: P = 3.51 GPa
步骤 5000: P = 11.57 GPa
步骤 6000: P = 14.89 GPa
步骤 7000: P = 5.96 GPa
步骤 8000: P = 6.58 GPa
高压平衡态结果:
最终压力: 10.06 GPa
体积变化: 8.96%
密度变化: 7.52%
估算体积模量: 112.3 GPa
# 高压模拟结果可视化
fig = make_subplots(
rows=2,
cols=2,
subplot_titles=("压力收敛", "体积压缩", "密度增加", "P-V关系"),
vertical_spacing=0.15,
horizontal_spacing=0.12,
)
# 压力演化
fig.add_trace(
go.Scatter(
x=hp_data["time"],
y=hp_data["pressure"],
mode="lines",
name="P(t)",
line=dict(color="blue", width=1),
),
row=1,
col=1,
)
fig.add_hline(
y=10,
line_dash="dash",
line_color="red",
annotation_text="目标: 10 GPa",
row=1,
col=1,
)
# 体积演化
fig.add_trace(
go.Scatter(
x=hp_data["time"],
y=hp_data["volume"],
mode="lines",
name="V(t)",
line=dict(color="green", width=1),
),
row=1,
col=2,
)
# 密度演化
fig.add_trace(
go.Scatter(
x=hp_data["time"],
y=hp_data["density"],
mode="lines",
name="ρ(t)",
line=dict(color="purple", width=1),
),
row=2,
col=1,
)
# P-V关系（状态方程）
fig.add_trace(
go.Scatter(
x=hp_data["volume"][::10],
y=hp_data["pressure"][::10],
mode="markers",
name="P-V",
marker=dict(size=3, color="orange"),
),
row=2,
col=2,
)
# 更新坐标轴
fig.update_xaxes(title_text="时间 (fs)", row=1, col=1)
fig.update_xaxes(title_text="时间 (fs)", row=1, col=2)
fig.update_xaxes(title_text="时间 (fs)", row=2, col=1)
fig.update_xaxes(title_text="体积 (ų)", row=2, col=2)
fig.update_yaxes(title_text="P (GPa)", row=1, col=1)
fig.update_yaxes(title_text="V (ų)", row=1, col=2)
fig.update_yaxes(title_text="ρ (amu/ų)", row=2, col=1)
fig.update_yaxes(title_text="P (GPa)", row=2, col=2)
fig.update_layout(title="高压NPT模拟结果 (300K, 10GPa)", height=700, showlegend=False)
pshow(fig)
4.7. 第五部分：总结与最佳实践
4.7.1. 5.1 MTK算法优势
相比其他压力控制方法，MTK算法具有以下优势：
严格的统计力学基础：产生正确的NPT系综分布
时间可逆性：对称Trotter分解保证长时间稳定性
正确的涨落：体积涨落与压缩率正确关联
灵活性：可处理各向异性压力控制
4.7.2. 5.2 参数选择指南
系统类型 |
dt (fs) |
τ_T (fs) |
τ_P (fs) |
链长 |
|---|---|---|---|---|
金属 |
0.1-0.5 |
20-100 |
500-1000 |
3 |
共价晶体 |
0.1-0.3 |
50-200 |
1000-2000 |
3-5 |
分子晶体 |
0.5-1.0 |
100-500 |
2000-5000 |
3-5 |
液体 |
1.0-2.0 |
100-500 |
1000-5000 |
3 |
4.7.3. 5.3 常见问题与解决方案
压力振荡过大
增大τ_P（压力弛豫时间）
增加系统尺寸
检查力计算精度
体积不收敛
延长平衡时间
调整初始体积更接近目标
检查势函数截断半径
守恒量漂移
减小时间步长
增加链长度
检查数值精度设置
4.7.4. 5.4 与其他方法比较
方法 |
优点 |
缺点 |
适用场景 |
|---|---|---|---|
Berendsen |
简单、快速收敛 |
不产生正确涨落 |
初始弛豫 |
Parrinello-Rahman |
可处理各向异性 |
可能不稳定 |
晶体相变 |
MTK |
严格、稳定 |
计算成本高 |
精确模拟 |
Monte Carlo |
严格采样 |
无动力学信息 |
平衡态性质 |
4.7.5. 5.5 进一步学习
原始文献：Martyna et al., J. Chem. Phys. 101, 4177 (1994)
各向异性压力控制
非平衡NPT模拟
与自由能计算结合
4.8. 练习与思考
4.8.1. 练习1：参数影响研究
修改上述代码，系统研究：
时间步长对守恒量的影响
链长度对温度/压力涨落的影响
系统尺寸对压力涨落的影响
4.8.2. 练习2：热膨胀系数计算
在不同温度下运行NPT模拟，计算铝的热膨胀系数：
\[\alpha = \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_P\]
4.8.3. 练习3：状态方程拟合
在不同压力下运行NPT模拟，拟合Birch-Murnaghan状态方程：
\[P = \frac{3B_0}{2}\left[\left(\frac{V_0}{V}\right)^{7/3} - \left(\frac{V_0}{V}\right)^{5/3}\right]\]
4.8.4. 思考题
为什么MTK算法需要压浴链？仅用热浴链是否足够？
如何判断NPT模拟是否达到平衡？
在什么情况下需要使用各向异性压力控制？
4.9. API参考
本教程中使用的主要类和函数：
4.9.1. 核心类
4.9.2. MD方案
MTKNPTScheme: MTK NPT积分方案create_mtk_npt_scheme(): MTK方案工厂函数
4.9.3. 传播子
MTKBarostatPropagator: MTK压力控制器NoseHooverChainPropagator: Nosé-Hoover链恒温器
4.9.4. 势函数
EAMAl1Potential: 铝EAM势
4.9.5. 力学计算
StressCalculator: 应力张量计算
详细文档请参考 API文档。