LAMMPS计算表面张力：从力学定义到多方法对比

表面张力是液体界面的核心热力学性质，在分子动力学模拟中可通过多种路径获得。LAMMPS提供了从压力张量分析到界面构型拟合等不同层次的方法，每种方法的适用场景和精度边界各不相同。

表面张力的力学定义

在分子模拟中，表面张力γ的定义来自Kirkwood-Buff公式：

γ = (1/2A) * ⟨Σ_ij (r_ij – 3*r_ij,z²/r_ij) * φ'(r_ij)⟩

其中A是界面面积，求和遍及所有分子对，r_ij,z是分子间距离的z分量，φ'(r_ij)是势能导数。这个定义的本质是：表面张力等于界面区域法向压力与切向压力之差的时间平均。

LAMMPS中实际计算的是机械压力张量，表面张力通过各向异性压力分量获得：

γ = L_z * (P_zz – (P_xx + P_yy)/2)

L_z是模拟盒子z方向长度，P_xx、P_yy是切向压力（平行于界面），P_zz是法向压力（垂直于界面）。对于平板型液-气界面模型，P_zz接近零（平衡时法向压力≈0），因此γ ≈ -L_z * (P_xx + P_yy)/2。

方法一：压力张量法（最直接）

这是最经典也最常用的方法。构建一个沿z方向延伸的液-气界面平板模型，LAMMPS自动计算压力张量的三个对角分量。

模型构建

以水分子（SPC/E力场）为例：

• 在x-y平面构建矩形液相区域（如30Å×30Å）
• z方向留出足够真空（总z长度60-80Å），使液相两侧各形成气-液界面
• 周期性边界条件：x/y方向周期，z方向也周期（形成上下两个界面）

LAMMPS脚本核心

# 力场设置
pair_style lj/cut/coul/long 10.0 10.0
pair_coeff 1 1 0.1553 3.166 10.0
bond_style harmonic
angle_style harmonic
kspace_style pppm 1.0e-4

# 界面建立后平衡
fix 1 all nvt temp 298.0 298.0 100.0
run 200000  # 1ns平衡

# 切换到NVE+恒温器（NVT会干扰压力张量统计）
unfix 1
fix 1 all nve
fix 2 all langevin 298.0 298.0 100.0 12345
run 1000000  # 5ns生产

表面张力提取

LAMMPS运行中输出压力分量：

thermo_style custom step temp press pxx pyy pzz
thermo 1000

表面张力的后处理公式：

γ = (L_z/2) * <(P_xx + P_yy)/2 – P_zz>

除以2是因为有两个界面（上下各一个）。注意L_z是平衡后的实际盒子长度。

SPCE水的实测值

用上述方法，SPC/E水在298K下的表面张力典型值为：

• 2ns平衡 + 5ns生产：γ = 62.8 ± 2.1 mN/m
• 实验值：71.97 mN/m
• 低估约13%，这是SPC/E力场的已知偏差（SPC/E的修正参数针对体相性质优化，界面性质有系统性偏低）

方法二：测试面积法（TAM）

当压力张量法的统计噪声太大（如界面面积小或密度低）时，测试面积法提供更稳定的替代方案。

原理是对平衡后的构型做各向同性面积缩放（保持分子数和体积不变），计算缩放前后的能量差：

γ = (1/2) * (dE/dA) * A

LAMMPS实现：

# 用fix deform做面积扫描
variable s equal 0.98
fix stretch all deform 1 x scale ${s} y scale ${s} remap x
run 0  # 不做动力学，只算能量

对一系列缩放因子（0.95-1.05）重复计算，拟合E vs A曲线的斜率即为表面张力。这种方法的优势是对力场中长程校正更敏感，劣势是需要在多个面积点重复计算。

方法三：Young-Laplace方程法（适用于液滴）

对于球形液滴模型，可以用Young-Laplace方程反推表面张力：

ΔP = 2γ/R

在LAMMPS中构建球形液滴，统计液滴内部和外部的压力差ΔP，以及液滴半径R，即可计算γ。

# 构建球形液滴
region drop sphere 0 0 0 25
create_atoms 1 region drop

# 平衡后统计径向压力分布
compute rdf all rdf 100
fix out all ave/time 100 1 100 c_rdf[*] file pressure_profile.txt mode vector

这种方法的适用场景是纳米尺度液滴（R<10nm），此时曲率效应显著，平板模型不再适用。但需要注意Gibbs-Tolman修正：纳米液滴的表面张力比平面界面偏低约10-20%。

力场选择对结果的影响

不同力场对水的表面张力计算结果差异显著：

TIP4P/2005是目前纯水界面模拟的首选力场，精度和计算效率平衡最佳。如果体系涉及离子或有机分子，需要验证该力场在混合体系下的表现。

统计精度的关键因素

表面张力是涨落量，统计误差与模拟时间的关系为σ ~ 1/√(t)。实际经验：

• 液-气界面：需要至少2-5ns生产时间，统计误差控制在±2 mN/m
• 液-液界面：需要5-10ns，因为界面波动更慢
• 界面面积越大统计越稳定，但计算量也线性增长

一个常见错误是平衡时间不够——液相密度未收敛到平衡值时就开始生产统计，会导致表面张力偏移5-10 mN/m。判断平衡是否充分的标志是液相体密度在100ps以上保持稳定波动。