LAMMPS计算自由能 — 从热力学积分到伞形采样的实战方法

LAMMPS计算自由能与GROMACS有所不同——LAMMPS更侧重于材料科学体系（金属、陶瓷、高分子），而非生物分子。LAMMPS计算自由能的典型场景包括：空位形成自由能、晶界迁移势垒、离子在固体电解质中的迁移自由能、液体-固体界面吸附自由能等。LAMMPS本身不内置完整的自由能计算框架，需要借助COLVARS模块、PLUMED插件或自定义脚本实现。

一、LAMMPS自由能计算方法概览

LAMMPS支持的自由能计算方法：

伞形采样（COLVARS）：COLVARS是LAMMPS内置的集合变量模块，可以定义反应坐标（距离、角度、RMSD等）并施加谐调约束。适合一维PMF计算。

热力学积分（TI）：通过逐步改变势能参数（λ从0到1）计算自由能差。LAMMPS的fix adapt命令可以动态修改势能参数，实现alchemical变换。

 Metadynamics：通过PLUMED插件实现。在已有轨迹上不断添加高斯势，迫使体系离开当前势阱探索新的构型空间。适合多维自由能面计算。

Einstein关系/Green-Kubo：不直接算自由能，而是通过扩散系数（Einstein）或涨落积分（Green-Kubo）间接获取热力学量。适合计算输运性质。

二、COLVARS伞形采样配置

COLVARS是LAMMPS中最方便的自由能计算工具。配置方式：

在LAMMPS输入脚本中添加：

# 加载COLVARS模块
fix 1 all colvars colvars_config.in

# 输出COLVARS数据
fix 1 all colvars colvars_config.in output colvars_output.txt

COLVARS配置文件（colvars_config.in）：

colvar {
name dist
distance {
group1 { atomNumbers 1-50 }
group2 { atomNumbers 51-100 }
}
lowerBoundary 5.0
upperBoundary 20.0
width 0.2    # 窗口宽度
}

harmonic {
colvars dist
forceConstant 5.0    # kcal/mol/Ų
centers 8.0    # 目标距离
}

每个窗口需要单独运行一次模拟，修改centers值。或者使用COLVARS的窗口扫描功能自动跑多个窗口。

三、热力学积分方法

LAMMPS的热力学积分通过fix adapt实现。以计算空位形成自由能为例：

目标：计算从完整晶体中移除一个原子（A→vacuum）的自由能变化。

策略：不直接删除原子，而是逐步减小该原子的LJ相互作用参数，使其”幽灵化”。

# fix adapt逐步改变原子1的epsilon
fix 1 all adapt 1 pair lj/cut epsilon 1 1 v_lambda

其中v_lambda是一个变量，从1.0线性变到0.0。需要设置20-40个λ值，每个值做独立模拟。

每个λ值下计算⟨∂U/∂λ⟩，然后对λ从0到1积分：ΔG = ∫⟨∂U/∂λ⟩dλ。

LAMMPS不自动输出∂U/∂λ，需要用fix print输出势能分量，后处理计算。或者用PLUMED的THERMODYNAMIC_INTEGRATION模块自动处理。

四、实战复盘：Li离子在Li₃PS₄固体电解质中的迁移自由能

体系：Li₃PS₄固体电解质，正交相，超胞3×3×2（共648原子）。目标计算Li离子从一个晶格位迁移到相邻空位的迁移自由能曲线（PMF）。

反应坐标：Li离子与空位之间的距离。初始位置3.2 Å（晶格位），终点0 Å（到达空位位置，即离子跨过势垒到达新位点）。

第一步：用COLVARS定义反应坐标。

colvar {
name li_vac_dist
distance {
group1 { atomNumbers 145 }    # 迁移Li离子
group2 { dummyAtom 0 0 0 }    # 空位中心坐标
}
width 0.1
}

第二步：窗口设置。沿迁移路径设20个窗口，间距0.15 Å，力常数3.0 kcal/mol/Ų。每个窗口在LAMMPS中用NVT系综跑2 ns（500 K，加速扩散），前0.5 ns平衡，后1.5 ns采样。

第三步：WHAM重构。用COLVARS自带的WHAM工具或外部WHAM程序处理20个窗口的力数据，得到PMF曲线。

第四步：结果分析。PMF显示一个明显的势垒，峰高0.52 eV（50.1 kJ/mol），位于迁移路径中点（距离1.6 Å处）。势垒值与实验测量的活化能0.49 eV吻合（偏差6%）。

关键发现：迁移路径不是直线，而是经过一个中间位（tetrahedral位点）。PMF在中间位有一个浅谷（深度0.08 eV），意味着Li离子在迁移过程中会短暂停留。这个两步迁移机制与第一性分子动力学的结果一致。

第五步：温度依赖分析。在400 K、450 K、500 K、550 K四个温度下重复计算，提取势垒高度的温度依赖性。结果：势垒高度随温度升高略有下降（从0.55 eV降到0.49 eV），表明迁移是热激活过程，活化能随温度有约0.06 eV的变化。

五、常见问题与解决方法

COLVARS窗口重叠不足：如果相邻窗口的直方图没有重叠区域，WHAM重构会出现缝隙，PMF曲线不连续。检查方法：绘制所有窗口的直方图叠加图，确认相邻窗口在峰值一半高度处有重叠。如果不重叠，增大力常数或减小窗口间距。

COLVARS力常数单位混淆：LAMMPS的默认单位是real（kcal/mol/Ų），但很多教程用metal单位（eV/Ų）。如果单位设置错误，力常数可能差24倍（1 eV = 23.06 kcal/mol）。务必在COLVARS配置中检查units设置。

LAMMPS输出与WHAM的格式对接：COLVARS输出的力数据格式可能与某些WHAM程序不兼容。最安全的方式是用COLVARS内置的runAve输出格式，或手动转换为GROMACS的pullf.xvg格式后用gmx wham处理。

固相采样的特殊性：固体中原子振动幅度小，采样效率低于液体。同一个窗口可能需要更长的平衡和产出时间（固体10-20 ns vs 液体2-5 ns）。判断标准：看COLVARS变量的时间序列是否在窗口中心附近做稳定振荡，如果还在漂移说明未平衡。

六、专业计算服务推荐

需要LAMMPS自由能计算服务？

科研学术网提供专业的LAMMPS自由能计算服务，涵盖COLVARS伞形采样、热力学积分、PLUMED Metadynamics、PMF重构全流程。

• ✅ 博士级工程师团队，一对一技术支持
• ✅ 计算结果可靠，可提供详细的技术报告
• ✅ 周期灵活，加急项目最快3天交付
• ✅ 价格透明，无隐形费用

立即咨询报价 →