LAMMPS计算自由能：自由能微扰与积分方法

LAMMPS计算自由能概述

自由能是热力学中描述体系稳定性和相平衡的核心物理量。在分子模拟中，自由能计算用于预测材料相稳定性、结合亲和力、溶解自由能和反应自由能等关键性质。LAMMPS作为通用分子动力学软件，通过多种方法可以实现自由能的计算，包括自由能微扰（FEP）、热力学积分（TI）和伞形采样等。

自由能计算在药物设计、材料筛选和相平衡研究中具有不可替代的地位，是连接微观模拟与宏观热力学的桥梁。

自由能计算的理论基础

自由能类型

自由能差计算

自由能差是自由能计算的核心目标：

ΔF = F(λ=1) – F(λ=0)

其中λ为耦合参数，控制从参考态(λ=0)到目标态(λ=1)的转变。

LAMMPS自由能计算方法

方法一：自由能微扰（FEP）

基本原理： ΔF = -kT·ln⟨exp[-β·ΔU]⟩_λ

其中ΔU为微扰前后的势能差，⟨…⟩_λ表示在λ状态下的系综平均。

LAMMPS实现：

FEP方法需要在多个λ值下进行独立模拟：

# LAMMPS FEP计算框架
# 对每个λ值：

# 1. 运行λ状态下的平衡MD
fix 1 all nvt temp 300 300 100

# 2. 采集λ'状态（微扰态）的能量
# 使用compute命令计算微扰势能
compute fep all fep 300 pair 1 1 v_lambda_s 1 0.0
# 计算从当前λ到λ'=0的微扰

# 3. 运行生产阶段
run 1000000

# 4. 后处理：BAR/FEP分析

FEP注意事项：

• 相邻λ值之间能量差应 < 2kT
• λ窗口数量：10-20个
• 每个窗口需要充分平衡
• 用BAR（Bennett Acceptance Ratio）方法提高精度

方法二：热力学积分（TI）

基本原理： ΔF = ∫₀¹ ⟨∂U/∂λ⟩_λ dλ

通过计算每个λ值下∂U/∂λ的系综平均，然后数值积分得到自由能差。

LAMMPS实现：

# LAMMPS TI计算框架
# 定义λ相关的势能
variable lambda equal 0.5  # 当前λ值

# 修改势能使其依赖λ
# U(λ) = (1-λ)·U_A + λ·U_B
# ∂U/∂λ = U_B - U_A

# 计算∂U/∂λ
compute dU all pair 1 1  # 计算势能导数

# 在NVT系综下运行
fix 1 all nvt temp 300 300 100
run 1000000  # 生产阶段

# 输出∂U/∂λ的时间序列
fix 2 all ave/time 100 10 1000 c_dU file dU_lambda_0.5.dat

多窗口TI流程：

1. 选择λ值：λ = 0, 0.1, 0.2, …, 0.9, 1.0
2. 对每个λ值运行独立MD
3. 提取每个λ的⟨∂U/∂λ⟩
4. 数值积分：梯形法或Simpson法
5. ΔF = ∫⟨∂U/∂λ⟩dλ

方法三：伞形采样（Umbrella Sampling）

基本原理： 通过施加偏置势沿反应坐标采样，然后用WHAM方法重构自由能面。

LAMMPS实现：

# 伞形采样设置
# 反应坐标：两粒子间距离r

# 定义反应坐标
compute 1 all pair/local dist
variable r equal c_1

# 施加偏置势（谐振势）
# U_bias = 0.5*k*(r - r₀)²
fix 1 all spring couple 1 2 10.0 0.0 0.0 0.0 5.0 0.0 0.0 0.0
# 参数：k=10.0, r₀=5.0

# 运行采样
fix 2 all nvt temp 300 300 100
run 1000000

# 不同窗口使用不同r₀值
# r₀ = 2, 3, 4, 5, 6, 7, 8 Å

WHAM后处理：

# 使用WHAM程序处理多窗口数据
# 输入：每个窗口的r直方图
# 输出：自由能面 F(r)
wham r_min r_max num_bins tol temperature metadata.dat output.dat

方法四：非平衡动力学（Jarzynski等式）

基本原理： ΔF = -kT·ln⟨exp(-β·W)⟩

其中W为非平衡过程所做的功。

LAMMPS实现：

# 快速拉伸/压缩模拟
fix 1 all nvt temp 300 300 100

# 施加随时间变化的力
variable f equal 10.0*dt
fix 2 all addforce 0 0 v_f

# 记录做功
# W = ∫F·dr
# 需要多个独立轨迹

LAMMPS自由能计算插件

FEP/TI插件（USER-FEP包）

LAMMPS的USER-FEP包提供了专门的自由能计算工具：

# 编译USER-FEP
cd lammps/src
make yes-USER-FEP
make mpi

compute fep命令：

compute fep all fep temp pair 1 1 v_lambda_s 1 0.0
# 参数：
# temp：温度
# pair：势能对类型
# v_lambda_s：λ参数变量
# 1 0.0：微扰到λ'=0.0

fix adapt/fep命令：

# 随模拟时间改变λ
fix 1 all adapt/fep 1 pair 1 1 v_lambda_s 0.0 1.0
# 从λ=0.0渐变到λ=1.0

PLUMED插件

PLUMED是强大的自由能计算增强采样插件，与LAMMPS兼容：

LAMMPS中调用PLUMED：

# 在LAMMPS输入文件中
fix 1 all plumed plumedfile input.dat

PLUMED输入文件（input.dat）：

# 定义反应坐标
d: DISTANCE ATOMS=1,2

# 伞形采样
u: BIASVALUE ARG=d

# 元动力学(Metadynamics)
METAD ARG=d SIGMA=0.5 HEIGHT=1.0 PACE=500

# 输出
PRINT ARG=d,u STRIDE=100 FILE=colvar.dat

PLUMED支持的方法：

• 伞形采样 + WHAM
• 元动力学（Metadynamics）
• 加速采样（ABMD）
• 变分增强采样（VES）

不同体系的自由能计算

1. 结合自由能

分子/颗粒结合自由能：

LAMMPS FEP方案：

# 阶段1：关闭静电（λ_chrg: 1→0）
# 阶段2：关闭LJ（λ_lj: 1→0）
# ΔG_bind = ΔG_chrg + ΔG_lj

# 每个阶段10个λ窗口

2. 相变自由能

固-液相变：

3. 表面自由能

晶体表面自由能： γ = (E_slab – N·E_bulk) / (2A)

# LAMMPS计算
# 1. 优化体相结构 → E_bulk
# 2. 构建slab → 优化 → E_slab
# 3. 计算表面积A
# 4. γ = (E_slab - n·E_bulk)/(2A)

4. 缺陷形成自由能

空位形成自由能： ΔF_vac = F(defect) – F(perfect) + μ_atom

# TI方法计算空位形成能
# λ=0：完美晶体
# λ=1：含空位晶体
# ΔF = ∫⟨∂U/∂λ⟩dλ

计算精度与误差控制

统计误差估算

收敛性检查

1. 时间收敛：检查ΔF随模拟时间的变化
2. 窗口收敛：检查每个λ窗口是否充分采样
3. hysteresis检查：正向(0→1)和反向(1→0)结果对比
4. block average：分块平均估算标准误差

常用误差标准

常见问题与解决方案

问题1：FEP/TI不收敛

原因：λ窗口间距过大或采样不足。

解决方案：

• 增加λ窗口数量（从10增加到20）
• 延长每个窗口的模拟时间
• 使用软核势避免采样困难
• 检查相空间重叠

问题2：正反向结果不一致（hysteresis）

原因：体系存在滞后效应或采样不完全。

解决方案：

• 使用BAR方法（自动处理正反向）
• 增加模拟时间
• 使用更密的λ窗口
• 检查是否存在相变

问题3：伞形采样窗口重叠不足

原因：相邻窗口的采样区域不重叠。

解决方案：

• 减小窗口间距
• 增加偏置势刚度
• 添加中间窗口
• 检查直方图重叠

问题4：自由能面不光滑

• 增加采样时间
• 增加窗口数
• 使用更小的bin宽度
• 检查反应坐标选择是否合理

实操案例：离子水合自由能

1. 体系：Na⁺在水中（1个Na⁺ + 999个水分子）
2. 方法：TI，分两阶段
   • 阶段1：关闭Na⁺电荷（λ_chrg: 1→0），10个窗口
   • 阶段2：关闭Na⁺-水LJ（λ_lj: 1→0），10个窗口
3. 设置：
   • NVT, T=298K, 每窗口5ns
   • 力场：SPC/E水 + Joung-Cheatham离子
4. 结果：
   • ΔG_chrg = -82.3 kcal/mol（电荷贡献）
   • ΔG_lj = -2.1 kcal/mol（LJ贡献）
   • 总水合自由能 = -84.4 kcal/mol
   • 实验值：-87.1 kcal/mol
   • 偏差：3.1%（合理范围）
5. 误差分析：
   • 标准误差：±0.5 kcal/mol
   • Hysteresis：< 0.3 kcal/mol

实操案例：颗粒聚集自由能

1. 体系：两个纳米颗粒在溶液中的聚集
2. 方法：伞形采样 + WHAM
3. 反应坐标：颗粒间距r（5-30 Å）
4. 窗口：20个，间距1.25 Å
5. 偏置势：k=5.0 kcal/mol/Ų
6. 设置：NVT, T=300K, 每窗口2ns
7. 结果：
   • 自由能极小值在r=12 Å
   • 聚集自由能ΔG = -15.2 kcal/mol
   • 能垒：3.5 kcal/mol（解离需要克服）

总结

LAMMPS通过FEP、TI、伞形采样和元动力学等方法，可以精确计算各类体系的自由能差。在实际应用中，需要根据体系特点选择合适的方法，合理设置耦合参数窗口，并通过收敛性检查和误差分析确保结果可靠性。结合PLUMED等增强采样插件，可以进一步扩展自由能计算的能力和应用范围。

我们提供专业的LAMMPS自由能计算服务，涵盖结合自由能、溶解自由能、相变自由能和缺陷形成自由能等各类计算需求，支持从力场设置到结果分析的完整计算流程。