GROMACS 二维密度与三维空间分布

同一条轨迹可以回答不同问题：沿膜法向哪里有水、膜中间哪些xy位置有水、水在三维空间停留在哪里，以及脂质尾链的取向怎样改变。本篇使用膜孔道建模教程中的真实Martini3 DPPC轨迹，把这些问题分别处理。

下载完整轨迹和分析脚本，解压后先执行：

bash view.sh

analysis/figures/pore-native-density-spatial.png 是下方的三联图。完整生产轨迹201帧、相邻5 ps，总长1 ns；W是粗粒化水珠。

1 选择什么数据和原子组#

参考目录有TPR、XTC、EDR、analysis-index.ndx、cylinder-definition.json。索引中的 Water 为全部W珠，Bilayer 为全部DPPC，Phosphates 为PO4，TailVectors 按每条脂质C1A、C4A顺序排列。

cd reference-results/v100-pore-r2p4-s13579
head -n 16 analysis-index.ndx
cat cylinder-definition.json
module load gromacs/2026.4-v100-threadmpi
mkdir manual-analysis

以下命令输出到新建的 manual-analysis，保留原始参考输出。已有此目录时换一个新名称。

2 沿z方向的数密度：gmx density#

printf 'Bilayer
Water
' | gmx density -s md.tpr -f md.xtc   -n analysis-index.ndx -d Z -sl 80 -dens number -center   -o manual-analysis/water-z-density.xvg

先选用于居中的膜，再选要统计的水。-sl 80 将法向分为80层；-dens number 输出W珠数密度，单位nm⁻³。-center 使图中坐标相对膜中心，不直接当作未平移的盒子绝对z坐标。膜中间密度低说明这个区域少水，不等于测得自由能障碍或水通透率。

3 固定膜中平面薄层的二维图：gmx densmap#

先从实际圆柱定义获取中平面上下边界：

python3 - <<'PY'
import json
d=json.load(open('cylinder-definition.json'))
print('zmin =',d['center_nm'][2]-d['halfheight_nm'])
print('zmax =',d['center_nm'][2]+d['halfheight_nm'])
PY

把输出值填到下方的 ZMIN、ZMAX。可用下面两行直接读取同样的值：

ZMIN=$(python3 -c 'import json;d=json.load(open("cylinder-definition.json"));print(d["center_nm"][2]-d["halfheight_nm"])')
ZMAX=$(python3 -c 'import json;d=json.load(open("cylinder-definition.json"));print(d["center_nm"][2]+d["halfheight_nm"])')
printf 'Water
' | gmx densmap -s md.tpr -f md.xtc -n analysis-index.ndx   -aver z -n1 32 -n2 32 -unit nm-2 -xmin "$ZMIN" -xmax "$ZMAX"   -od manual-analysis/midplane-water-nm2.dat   -o manual-analysis/midplane-water.xpm -dmax 4

-aver z 沿z方向投影，-xmin/-xmax 在这里选择投影方向的薄层边界，并不是限制图上的x轴。图上每个横向网格统计该薄层内的水。

本例选 nm-2，数据是平均W珠面积密度。网格32×32，要得到每格平均W珠数，乘以 Lx*Ly/(32*32)。不要把nm⁻²误写成nm⁻³，也不要把本版本 -unit count 输出不经核对就当作每格平均粒子数。

.dat 第一行/列是坐标，内部32×32为数值。XPM适合快速查看，定量处理使用DAT。附件脚本直接读取数值重新绘图。

4 三维占据场：gmx spatial#

cd manual-analysis
printf 'Water
Bilayer
' | gmx spatial -s ../md.tpr -f ../md.xtc   -n ../analysis-index.ndx -bin .25 -nodiv -nab 8
cd ..

第一组选要统计的水，第二组选同时写入CUBE结构区的膜。grid.cube 生成在当前工作目录，所以单独进入输出目录。-bin .25 为0.25 nm体素；-nodiv 得到每帧每体素的平均W珠数，不是归一到体相的分布函数，也不是原子级电子密度。

CUBE头部的坐标和网格向量使用Bohr，换算为nm时乘0.05291772109。不能直接把头部的4.724315看成4.7 nm，它对应约0.25 nm。

本例不对蛋白质拟合，使用固定盒子的空间网格。用于溶剂化壳层时，应先解决分子完整性、PBC跳跃及拟合参考，否则溶质平移会抹平空间分布。这里没有蛋白质，也不把膜水图称为蛋白溶剂化壳。

当前原生输出的网格范围省略了每个方向最外一层占据体素；边缘部分水珠因此不在导出的CUBE内。附件分析显式计算这个范围与遗漏占据数，不把CUBE总和强行解释为体系全部水珠数。

5 尾链取向：gmx gangle与P2#

gmx gangle -s md.tpr -f md.xtc -n analysis-index.ndx   -g1 vector -g2 z -group1 'group "TailVectors"'   -oall manual-analysis/tail-angles.xvg   -oav manual-analysis/tail-angle-mean.xvg -xvg none

本例定义从C1A指向C4A的粗粒化向量，与膜法向z比较夹角θ，计算 P2=(3*cos(theta)^2-1)/2。接近平行法向时趋向1，垂直时为−0.5。跨周期边界的向量要按最小镜像处理。

这个P2不是原子级C–D氘序参量。 gmx order 面向适当的饱和联合原子链等输入，不直接用于这套Martini粗粒化珠。原子级尾链需按实际拓扑与氢原子定义另行准备，不能把珠子改名就当作氘序参量。

6 固定区域的水珠数：gmx select#

先查看程序按完整精度TPR生成的选择表达式：

python3 - <<'PY'
import json
print(json.load(open('cylinder-definition.json'))['native_selection'])
PY
SELECT=$(python3 -c 'import json;print(json.load(open("cylinder-definition.json"))["native_selection"])')
gmx select -s md.tpr -f md.xtc -select "$SELECT"   -os manual-analysis/water-cylinder-count.xvg

表达式对应固定轴心、半径0.5 nm、总高1 nm的圆柱。时间序列有201个保存时刻。更换体系应重定轴心和尺寸；孔移动时，固定圆柱只能回答原位置的占据情况。

7 把这些图放在一起解释#

一维图看水沿法向的分层，二维图定位中平面水占据的横向位置，三维图显示空间分布，圆柱计数给出固定区域的时间变化，尾链P2描述取向。它们观察同一条短轨迹的不同方面，不能相互替代。

本例8组全部保留：较大初始缺口在生产开始尚有水，之后观察区域大多变干；较小缺口在生产起点已干。结果可用于理解建模、选区和后处理流程，不能据此宣称稳定开放孔、测得通透率或获得平衡材料参数。

参数说明见GROMACS densmap手册、spatial手册、gangle手册和order输入说明。