MDAnalysis水的几何、密度与介电响应

同一份水轨迹可以回答三个不同问题：O—H距离是否正确、氧原子在空间如何分布、总偶极矩如何涨落。本文逐步展示终端操作，并说明三个结果各自的单位和适用范围。

1 下载和运行#

下载水、泛素与珠簧聚合物的实际轨迹、分析脚本和完整数据，保存为mdanalysis-water-polymer-wave330-v2.zip：

module load pymatgen/2026.9.24
python -m zipfile -e mdanalysis-water-polymer-wave330-v2.zip .
cd mdanalysis-water-polymer
bash run.sh
source env.sh

run.sh只对已完成的真实轨迹做后处理，生成results/和六张plots/图片，不启动新的分子动力学。请在新的解压目录运行，避免覆盖自己的输出。MDAnalysis为2.9.0，环境通过现有qvasp 2.26的Python选择，没有调用qvasp或VASP；其他机器可通过ANALYSIS_PYTHON指定兼容环境运行run.sh。

水体系来自已有GROMACS水模拟，884个TIP3P水、2652原子、1001帧，时间0–100 ps。水分析每20帧取一帧，共51帧。轨迹为NPT，实际盒长约29.87494–30.29525 Å，不能当作固定体积轨迹交给只支持固定正交盒的LinearDensity。本例使用能明确处理所需空间网格或平均体积的工具；短轨迹统计仍要检查采样收敛。

2 AtomicDistances：逐对O—H距离#

python - <<'PYCODE'
import MDAnalysis as mda
from MDAnalysis.analysis.atomicdistances import AtomicDistances
u = mda.Universe('inputs/water/water.gro', 'inputs/water/water.xtc')
o = u.select_atoms('name OW')[:20]
h = u.select_atoms('name HW1')[:20]
a = AtomicDistances(o, h, pbc=True).run(stop=1001, step=20)
print(a.results.shape, a.results.mean())
PYCODE

本例得到(51,20)及平均 0.9574386 Å。第1个O对应第1个HW1，依次配对，不是20×20的所有距离；两个选择必须长度相同。pbc=True按每帧实际盒子处理最小镜像，水分子跨过边界时不能直接拿原始坐标相减。

3 DensityAnalysis：固定空间网格里的氧密度#

python - <<'PYCODE'
import numpy as np
import MDAnalysis as mda
from MDAnalysis.analysis.density import DensityAnalysis
u = mda.Universe('inputs/water/water.gro', 'inputs/water/water.xtc')
a = DensityAnalysis(u.select_atoms('name OW'), delta=1.0, padding=0,
                    gridcenter=[15.15]*3, xdim=36, ydim=36, zdim=36)
a.run(stop=1001, step=20)
d = a.results.density
print(d.grid.shape)
print('mean oxygen count:', d.grid.sum()*np.prod(d.delta))
d.export('my-oxygen-density.dx')
PYCODE

实际网格为 37×37×37，应以d.edges和d.delta为准，不要把输入尺寸直接当成最终格点数。氧数密度单位为Å⁻³；本例对完整包围网格积分得到 884，与每帧水分子数一致。它不是质量密度，也不是在所有体积变化条件下都能直接平均成一个热力学密度的通用方法。

图中边缘衰减同时受选定空间网格及NPT盒子变化影响，不应当作水形成了表面。若研究蛋白周围的水密度，先把蛋白对齐到统一参照，再定义固定网格，并让轨迹覆盖足够构象。

4 介电分析前：电荷和完整水分子从哪里来#

GRO没有可靠的力场电荷或键表；仅加载GRO和XTC不足以计算总偶极涨落。本包从同一原生TPR读取一个SOL分子类型、884个分子和SETTLE原子三元组：O、H、H。质量为16、1.008、1.008，电荷为−0.834、+0.417、+0.417 e，每个水分子电中性。

bash export-topology.sh
cat topology-export/water-topology.json

该脚本实际使用gmx dump -s inputs/water/water.tpr提取数值。刚性水的连接约束存在SETTLE里，本例的PDB导出没有CONECT，不能据此认为水没有连接关系。下面根据这个明确三元组添加两条O—H连接，用于把每个水分子完整展开，不靠几何猜键。这个转换器限定为本例单一SOL分子类型；其他力场、水模型、混合物或离子体系须改用相应的真实拓扑。

5 DielectricConstant：总偶极涨落#

python - <<'PYCODE'
import json
import numpy as np
import MDAnalysis as mda
from MDAnalysis.analysis.dielectric import DielectricConstant
u = mda.Universe('inputs/water/water.gro', 'inputs/water/water.xtc')
t = json.load(open('water-topology.json'))
u.add_TopologyAttr('charges', np.tile(t['charges'], t['molecules']))
u.add_TopologyAttr('masses', np.tile(t['masses'], t['molecules']))
u.add_bonds([(3*i,3*i+j) for i in range(t['molecules']) for j in (1,2)])
a = DielectricConstant(u.atoms, temperature=300, make_whole=True)
a.run(stop=1001, step=20)
print(a.results.M)
print(a.results.eps, a.results.eps_mean)
PYCODE

51帧的三个方向估计值约 86.48247、100.81877、54.52909，标量平均 80.61011。三个方向差别明显，说明有限采样的各向异性仍很大；不能把这个平均数直接当作已收敛介电常数，更不能据它宣称模型与实验一致。

总偶极为sum(q_i r_i)，单位eÅ。这个原生实现按样本平均体积归一化，并使用导体边界条件对应的涨落公式；对本例NPT短样本，应把它作为后处理操作演示，后续用更充分的采样和合适系综检查统计误差。体系包含自由离子或带电片段时，不要直接沿用这套无自由电荷的公式。官方介电分析说明列出了相应条件。

完整逐帧偶极、体积和三维网格在results/FULL-FIELDS.npz中；下一篇讲指定原子对与氢键存活。

补充：串行与双进程密度后处理#

下载并行密度与构象集合示例，保存为ensemble-map.zip。此处另取同一真实水轨迹的前20个完整分子、201个连续帧，先做整分子与NoJump，再在实验室坐标网格统计氧密度；这是选定20个水分子的密度示范，不是全体884个分子的体相密度或界面稳态密度。

python3 -m zipfile -e ensemble-map.zip .
cd mdanalysis-ensemble-map
bash run.sh
bash verify.sh
bash plot.sh

使用MDAnalysis 2.9.0的原生后端接口：

serial = DensityAnalysis(oxygen, **grid).run(backend='serial')
parallel = DensityAnalysis(oxygen, **grid).run(
    backend='multiprocessing', n_workers=2, n_parts=2
)

两次使用完全相同的100 Å见方网格、2 Å网格间距和201个真实坐标帧。所有125000个网格点逐项相同，并分别与独立的三维直方图吻合；密度乘网格体积求和为20个水氧。这里只证明同一输入下结果一致，没有用这个小例子宣称双进程一定更快。n_workers=2确实会使用两个CPU工作进程，自己的大轨迹请放在已分配资源内处理。