LAMMPS计算密度：从均匀体相到界面密度分布

密度是分子动力学模拟中最基础也是最容易出错的物理量之一。在均匀体系中，密度计算直截了当——总质量除以总体积；但在非均匀体系（界面、纳米孔、液滴）中，密度的空间分布才是真正有价值的分析目标。LAMMPS提供了多种密度计算工具，关键在于选择合适的统计方法和空间分辨率。

均匀体系密度的基本计算

对于体相液体或固体，密度可以直接从模拟盒子尺寸和原子质量推算：

variable rho equal mass(all)/(vol*1.66054)
thermo_style custom step temp v_rho

这里mass(all)返回所有原子总质量（原子质量单位），vol返回盒子体积（ų），除以1.66054是将g/mol·Å⁻³转换为g/cm³。

这种全局密度计算简单直接，但有两个前提：体系必须已达到平衡，且盒子体积在NPT系综下已收敛。在NVT系综中，密度是输入参数而非输出，通常需要在NPT下先确定平衡密度，再切换到NVT做生产模拟。

空间密度分布：分bin统计

界面体系的核心分析是密度沿某一方向的分布。LAMMPS用compute chunk和fix ave/chunk组合实现：

液-气界面密度剖面

以水-气界面为例，沿z方向统计氧原子密度：

# 分bin设置：z方向每0.5Å一个bin
compute chunks all chunk/atom bin/1d z lower 0.5 units box

# 统计每个bin内氧原子数密度
# group oxygen type 1（水氧原子类型为1）
compute chunks_ox oxygen chunk/atom bin/1d z lower 0.5 units box
fix density_profile oxygen ave/chunk 100 10 1000 chunks_ox &
    density/mass file rho_z.txt

输出文件rho_z.txt包含每行：bin编号、z坐标中心、质量密度（g/cm³）。

界面密度剖面的典型形态：液相区密度恒定（约1.0 g/cm³），过渡区密度从90%降到10%（界面层），气相区密度趋近零。界面厚度定义为密度从90%降至10%的z区间，SPC/E水在298K下的界面厚度约3.5-4.0 Å。

径向密度分布

对于球形对称体系（液滴、纳米颗粒周围的溶剂化层），需要径向密度分布：

compute rdf all rdf 100 1 1
fix gr all ave/time 100 10 1000 c_rdf[*] file rdf.txt mode vector

g(r)本身不是密度而是对均匀分布的归一化因子，要获得实际数密度n(r)：

n(r) = ρ_bulk * g(r)

其中ρ_bulk是体相数密度（atoms/ų）。

NPT系综下的密度收敛

NPT模拟中密度是动态量，收敛判断需要注意几个要点：

密度涨落与盒子尺寸

NPT中盒子体积自由涨落，密度标准差与体系大小相关：

• 1000个水分子：σ_ρ/ρ ≈ 2-3%
• 8000个水分子：σ_ρ/ρ ≈ 1%
• 27000个水分子：σ_ρ/ρ ≈ 0.5%

小体系的密度涨落是物理性的，不代表模拟不稳定。但如果密度平均值在长时间内持续漂移（>1 ns趋势），说明力场参数或恒温恒压器设置有问题。

恒压器选择

fix 1 all npt temp 298.0 298.0 100.0 iso 1.0 1.0 1000.0

Nosé-Hoover恒温恒压器（fix npt）是标准选择。阻尼参数tdamp建议100-200步（100-200 fs），pdamp建议500-2000步。pdamp太短会导致盒子剧烈振荡，太长则平衡慢。

对于力场给出明显偏离实验密度的体系（如TIP3P水在1 atm下平衡密度约0.98 g/cm³而非1.00），这是力场本身的系统偏差，不应该通过调节压力来”修正”。

二维密度图与等密度面

有时一维密度剖面不够，需要二维甚至三维密度分布。LAMMPS支持二维分bin：

compute chunks2d all chunk/atom bin/2d x lower 1.0 y lower 1.0 units box
fix 2d_map all ave/chunk 100 10 1000 chunks2d density/mass file rho_xy.txt

输出一个二维密度网格，可以用Origin、Python matplotlib或OVITO做等高线图。典型应用：

• 纳米通道内流体密度分布（x-z平面，揭示壁面吸附层）
• 液滴在表面上的接触角分析（x-z平面密度等高线）
• 蛋白质周围水化层密度（径向+角度二维分布）

密度计算中的常见陷阱

周期性边界条件的影响

在平板型界面模型中，z方向周期性意味着有两个界面（上下各一个）。统计密度时需要确认只取一个界面区域，否则两个界面的数据叠加导致密度剖面变形。

bin宽度选择

bin太窄（<0.5Å）会导致每个bin内原子数太少，统计噪声大；bin太宽（>2Å）会抹平界面细节。经验值：液态体系0.5-1.0Å，固态体系1.0-2.0Å。

质量密度 vs 数密度

LAMMPS的fix ave/chunk支持density/mass（质量密度）和density/number（数密度）。多组分体系中两者不等价——质量密度受重原子主导，数密度对轻原子更敏感。分析溶剂化层结构时通常用数密度（可以分辨不同原子类型的空间分布），分析宏观性质时用质量密度。

时间步数和输出频率

fix ave/chunk的参数Nevery Nrepeat Nfreq决定了统计方式。100 10 1000表示每100步采样一次，重复10次，每1000步输出一个平均值。如果Nfreq太小（如100），统计窗口太短，密度剖面噪声很大。建议Nfreq至少等于Nrepeat*Nevery，且总统计时间覆盖至少100 ps。