MDAnalysis质量与电荷密度剖面

本例使用集群已有的Martini 3 DPPC膜体系：198个DPPC、3050个粗粒化水珠，共5426个相互作用珠子。选用真实20 ps NVT阶段的41帧，保存间隔0.5 ps，盒子始终为80×80×100 Å。这里展示模型质量和电荷密度，不把粗粒化珠子的质量分布当作实验原子质量密度。

1 下载并运行#

下载固定盒子密度剖面完整示例，保存为linear-density.zip：

python3 -m zipfile -e linear-density.zip .
cd mdanalysis-linear-density
bash run.sh
bash verify.sh
bash plot.sh

环境为MDAnalysis 2.9.0。输入是inputs/nvt.gro、inputs/nvt.xtc和topology-metadata.json。质量、电荷与连接来自实际GROMACS 2026.4的TPR导出，没有猜测：珠子质量为54或72 u，DPPC的NC3带+1 e、PO4带−1 e，整个体系净电荷为零。

2 先核对固定盒子，再计算#

LinearDensity的归一化使用初始化时的盒子，适用于固定正交盒子。本例逐帧核对了全部41帧，不使用会变盒的NPT轨迹，也没有把原轨迹盒长改成常数。

u.add_TopologyAttr('masses', actual_masses)
u.add_TopologyAttr('charges', actual_charges)
u.add_bonds(actual_bonds)
u.trajectory.add_transformations(unwrap(u.atoms))
profile = LinearDensity(u.atoms, grouping='atoms', binsize=2.0).run()

binsize=2.0单位为Å。本例各盒长都能被2整除，切片体积为x/y方向16000 ų、z方向12800 ų。该版本用最长盒边建立共用坐标数组，所以三轴都输出50个区间；x/y中80 Å以后的条目是共用网格填充，不是盒中存在额外真空层。绘图和导出应按各轴真实长度截取。

3 原子、残基和连接片段有什么区别#

grouping='atoms'按每个粗粒化珠子的位置统计；residues按分子的质量中心统计；fragments按真实连接片段的质量中心统计。本例一个DPPC是一个残基和一个连接片段，水珠则各自单独成组，因此后两种结果相同。

先用真实连接把跨周期边界的分子恢复完整，再将相应分组包回盒子。仅把断开的分子坐标直接求平均，可能把质量中心放到错误位置。三种分组的空间剖面可以不同，但积分得到的总质量应相同。

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
f = np.load('results/FULL-FIELDS.npz')
for grouping in ['atoms', 'residues', 'fragments']:
    print(grouping, 'z mass density:', f[grouping + '_z_mass_density'][:5])
print('box:', f['actual_box_A'])
print('all model masses:', np.unique(f['model_masses_u']))
PYCODE

4 导出清楚的单位#

质量密度原生输出为g/cm³。为让电荷密度容易理解，下面按该版本的单位常数换成e/ų，同时导出质量密度：

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
from MDAnalysis.units import constants
f = np.load('results/FULL-FIELDS.npz')
e = f['edges_A']
center = (e[:-1] + e[1:]) / 2
rho_q = f['atoms_z_charge_density'] * constants['N_Avogadro'] * 1e-24
np.savetxt('z-density.dat', np.c_[center, f['atoms_z_mass_density'], rho_q],
           header='z_A mass_g_cm3 charge_e_A3')
print('total charge from z profile:', rho_q.sum() * f['slice_volumes_A3'][2])
PYCODE

下载包保留三轴质量/电荷的均值和标准差，以及各帧原始分箱权重。标准差描述这41帧的波动，不是均值的标准误差。残基/片段分组下每个DPPC净电荷为零，因此该分组的电荷剖面为零；这不代表分子内部的正负头基电荷消失了，珠子分组仍能显示它们。

这是20 ps热化阶段的使用示例，没有宣称膜结构或密度剖面已经平衡。自己的体系需要正确的质量、电荷、连接、固定盒子和足够长的采样，再解释密度层化与误差。

参考：MDAnalysis 2.9.0 LinearDensity接口。