4. 溶剂化与离子添加工程
完成拓扑构建后，体系仍是一个真空中的孤立分子。要让分子动力学模拟接近真实溶液环境，必须依次解决四个问题：用什么水、盒子多大、加多少离子、边界如何处理。这四步彼此耦合——水模型决定盒子尺寸与密度，盒子尺寸决定离子数量，而周期性边界条件（PBC）又反过来约束盒子的最小尺寸。本节按这条依赖链展开，并给出可直接套用的操作与判断依据。
一、水模型选择：不是越新越好，而是与力场匹配
水模型的核心差异在于其水分子几何与电荷分布被参数化的方式。常见三类：
- 三点刚性模型（SPC、SPC/E、TIP3P）：一个氧、两个氢、三个相互作用位点。计算量最小，与 AMBER、CHARMM、GROMOS 等生物力场兼容性最好。TIP3P 是 CHARMM36 的默认搭配，SPC/E 常用于 GROMOS 与部分通用力场。
- 四点模型（TIP4P、TIP4P/2005）：在氧原子沿 H-O-H 角平分线方向增加一个虚拟位点 M 承载负电荷，能更好再现水的密度与介电性质，但需要力场显式支持 M 位点，否则会与蛋白参数不兼容。
- 可极化与五点位模型（SWM4-NDP、TIP5P）：精度更高，但计算成本显著上升，且多数生物力场未配套参数。
工程判断原则：水模型必须与蛋白/配体力场同源或经文献验证兼容。若使用 CHARMM36 蛋白力场，选 TIP3P；若使用 GROMOS 54A7，选 SPC/E。混用（如 CHARMM36 + TIP4P/2005）会引入系统性误差，且虚拟位点需要额外的 virtual_sites 处理。本节示例统一采用 TIP3P + CHARMM36。
二、溶剂盒构建：形状、尺寸与最小镜像距离
盒子形状常用立方体（cubic）、菱形十二面体（dodecahedron）和截角八面体（truncated octahedron）。相同最小镜像距离下，菱形十二面体体积约为立方体的 0.71 倍，可显著减少水分子数，是球形蛋白的推荐选择。
尺寸确定的关键约束是最小镜像距离：溶质任意原子到盒子边界的距离，必须大于非键相互作用的截断半径（通常 1.0–1.2 nm）。工程上取溶质最大直径加两侧各 1.0–1.2 nm 缓冲。例如一个最大跨度约 4.0 nm 的蛋白，用菱形十二面体盒，盒边长可设为 4.0 + 2×1.2 ≈ 6.4 nm。
构建命令示例（以 GROMACS 为例）：
# 1. 定义盒子：菱形十二面体，溶质与盒边最小距离 1.2 nm
gmx editconf -f protein.gro -o box.gro \
-bt dodecahedron -d 1.2 -c
# 2. 填充 TIP3P 水
gmx solvate -cp box.gro -cs spc216.gro \
-o solv.gro -p topol.top执行后应检查输出：gmx solvate 会报告加入的水分子数，并自动在 topol.top 的 [ molecules ] 段追加 SOL 条目。若盒子过小，后续 gmx grompp 会报“minimum image distance too small”类错误，此时应回到 editconf 增大 -d。
三、离子中和与生理浓度添加
蛋白表面带电残基会使体系带净电荷。MD 中净电荷不为零会导致 PME 长程静电处理失效，因此必须先中和，再按目标浓度补加离子。
中和所需反离子数等于体系净电荷的绝对值。例如净电荷为 +8e，则需加入 8 个 Cl⁻。生理浓度通常取 0.15 mol/L（约 150 mM），对应在给定水体积中按 n = c × V × N_A 计算。GROMACS 的 genion 可一次性完成中和与浓度添加：
# 生成 tpr 以获取体系信息
gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
# 中和 + 补至 0.15 mol/L，替换 20 个水分子为离子
gmx genion -s ions.tpr -o ionized.gro -p topol.top \
-pname NA -nname CL -neutral -conc 0.15关键点：-neutral 先中和净电荷，-conc 0.15 再按 0.15 mol/L 补加 NaCl。genion 通过把溶剂分子替换为离子来实现，因此 topol.top 中 SOL 数量会相应减少，NA、CL 条目被追加。运行后应核对：净电荷是否为零、离子总数是否与目标浓度一致、水分子数是否合理减少。
gmx dump 或直接查看 topol.top 的 [ molecules ] 段，确认 NA 与 CL 数量之差等于体系净电荷的相反数；再用 gmx editconf -f ionized.gro -o check.pdb 目视检查离子是否均匀分布、有无重叠。
四、周期性边界条件设置
PBC 通过 .mdp 文件中的 pbc 参数控制，常用 xyz（三维周期）。与之配套的还有三组关键参数：
coulombtype = PME：长程静电用粒子网格 Ewald 方法，要求盒子最小镜像距离大于rlist。rvdw与rlist：范德华与邻居列表截断，通常取 1.0–1.2 nm，且rlist不小于rvdw。pbc = xyz配合nstlist邻居列表更新频率。
若模拟的是无限稀释的孤立分子（如小分子在真空或隐式溶剂中），可设 pbc = no，但此时不能使用 PME，需改用 cutoff 并配合足够大的盒子。对于溶液体系，始终使用 pbc = xyz + PME。
下面这张图把本节四个环节的依赖关系与关键参数串成一条工程链，便于对照检查：
五、常见错误与排查
实际操作中，以下三类问题最常出现：
| 现象 | 可能原因 | 处理方式 |
|---|---|---|
| grompp 报 “minimum image distance too small” | 盒子边长小于 2×rlist，或溶质贴近盒边 | 增大 editconf -d，或改用菱形十二面体 |
| genion 后体系仍带净电荷 | 未加 -neutral，或拓扑中残基电荷与预期不符 | 检查 -neutral，用 gmx dump 核对总电荷 |
| PME 报错或能量异常 | 盒子非中性，或 pbc = no 却用了 PME | 先中和；溶液体系保持 pbc = xyz |
-neutral 与 -conc 是协同工作的：先中和保证电中性，再按浓度补加。若手动分两步且顺序颠倒，容易在拓扑中留下净电荷，导致后续 PME 计算发散。
完成本节四步后，体系应满足：水模型与力场匹配、盒子最小镜像距离足够、净电荷为零、离子浓度符合目标、PBC 与 PME 配置一致。此时可进入能量最小化，检查是否存在原子重叠或异常高势能，再继续升温与平衡。