GROMACS 封端肽与力场实操
从PDB建拓扑，比较封端、端基质子化、AMBER/OPLS与水模型，实际提取phi/psi并画图。
阅读本篇目录
同一个名称“丙氨酸二肽”可能指不同的封端和末端状态。先明确化学组成、端基和电荷，再选力场、溶剂和离子参数，不能只靠改一个 -ff 选项就认为模型准备完成。
| 实际模型 | 力场 | 溶质原子 | 溶质电荷e |
|---|---|---|---|
| ACE-ALA-NME封端模型 | amber99sb-ildn | 22 | 0 |
| 非封端Ala-Ala | oplsaa | 23 | 0 |
| Ala-Ala的+1端基状态 | oplsaa | 24 | 1 |
| 单残基丙氨酸两性离子 | oplsaa | 13 | 0 |
封端模型初始几何来自 OpenMM8.2公开测试结构，随后由GROMACS按指定力场重新生成氢原子。Ala-Ala参考同一肽键几何另行构建，并经原生建模、最小化和平衡；配图取真实轨迹。它们不是相同化学组成，不能直接用总能量相减作为端基或质子化自由能。
1. 先生成溶质拓扑，再建立水溶液#
以封端ACE-ALA-NME为例：
gmx pdb2gmx -f alanine-dipeptide.pdb -o processed.gro \
-p topol.top -ff amber99sb-ildn -water tip3p -ignh
gmx editconf -f processed.gro -o box.gro -c -box 3 3 3 -bt cubic
gmx solvate -cp box.gro -cs spc216.gro -o solv.gro -p topol.top
-ignh让程序按所选库重新补氢。建模后先确认分子原子数和总电荷，再加离子与平衡。下面是本次原生生成的完整22原子GRO：
GROup of MAchos and Cynical Suckers
22
1ACE CH3 1 0.200 0.209 0.000
1ACE HH31 2 0.135 0.285 0.000
1ACE HH32 3 0.186 0.153 -0.082
1ACE HH33 4 0.186 0.153 0.082
1ACE C 5 0.343 0.264 -0.000
1ACE O 6 0.439 0.188 -0.000
2ALA N 7 0.356 0.397 -0.000
2ALA H 8 0.273 0.454 0.000
2ALA CA 9 0.485 0.461 -0.000
2ALA HA 10 0.528 0.429 0.084
2ALA CB 11 0.566 0.422 -0.123
2ALA HB1 12 0.655 0.468 -0.121
2ALA HB2 13 0.580 0.323 -0.124
2ALA HB3 14 0.517 0.450 -0.206
2ALA C 15 0.471 0.613 0.000
2ALA O 16 0.360 0.665 0.000
3NME N 17 0.585 0.683 0.000
3NME H 18 0.672 0.635 0.000
3NME CH3 19 0.585 0.828 0.000
3NME HH31 20 0.679 0.862 0.000
3NME HH32 21 0.537 0.862 -0.082
3NME HH33 22 0.537 0.862 0.082
0.54394 0.70872 0.28981
完整有效拓扑内容（仅省略自动生成的文件头注释）
; Include forcefield parameters
#include "amber99sb-ildn.ff/forcefield.itp"
[ moleculetype ]
; Name nrexcl
Protein 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
; residue 1 ACE rtp ACE q 0.0
1 CT 1 ACE CH3 1 -0.3662 12.01
2 HC 1 ACE HH31 2 0.1123 1.008
3 HC 1 ACE HH32 3 0.1123 1.008
4 HC 1 ACE HH33 4 0.1123 1.008
5 C 1 ACE C 5 0.5972 12.01
6 O 1 ACE O 6 -0.5679 16 ; qtot 0
; residue 2 ALA rtp ALA q 0.0
7 N 2 ALA N 7 -0.4157 14.01
8 H 2 ALA H 8 0.2719 1.008
9 CT 2 ALA CA 9 0.0337 12.01
10 H1 2 ALA HA 10 0.0823 1.008
11 CT 2 ALA CB 11 -0.1825 12.01
12 HC 2 ALA HB1 12 0.0603 1.008
13 HC 2 ALA HB2 13 0.0603 1.008
14 HC 2 ALA HB3 14 0.0603 1.008
15 C 2 ALA C 15 0.5973 12.01
16 O 2 ALA O 16 -0.5679 16 ; qtot 0
; residue 3 NME rtp NME q 0.0
17 N 3 NME N 17 -0.4157 14.01
18 H 3 NME H 18 0.2719 1.008
19 CT 3 NME CH3 19 -0.149 12.01
20 H1 3 NME HH31 20 0.0976 1.008
21 H1 3 NME HH32 21 0.0976 1.008
22 H1 3 NME HH33 22 0.0976 1.008 ; qtot 0
[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
1 3 1
1 4 1
1 5 1
5 6 1
5 7 1
7 8 1
7 9 1
9 10 1
9 11 1
9 15 1
11 12 1
11 13 1
11 14 1
15 16 1
15 17 1
17 18 1
17 19 1
19 20 1
19 21 1
19 22 1
[ pairs ]
; ai aj funct c0 c1 c2 c3
1 8 1
1 9 1
2 6 1
2 7 1
3 6 1
3 7 1
4 6 1
4 7 1
5 10 1
5 11 1
5 15 1
6 8 1
6 9 1
7 12 1
7 13 1
7 14 1
7 16 1
7 17 1
8 10 1
8 11 1
8 15 1
9 18 1
9 19 1
10 12 1
10 13 1
10 14 1
10 16 1
10 17 1
11 16 1
11 17 1
12 15 1
13 15 1
14 15 1
15 20 1
15 21 1
15 22 1
16 18 1
16 19 1
18 20 1
18 21 1
18 22 1
[ angles ]
; ai aj ak funct c0 c1 c2 c3
2 1 3 1
2 1 4 1
2 1 5 1
3 1 4 1
3 1 5 1
4 1 5 1
1 5 6 1
1 5 7 1
6 5 7 1
5 7 8 1
5 7 9 1
8 7 9 1
7 9 10 1
7 9 11 1
7 9 15 1
10 9 11 1
10 9 15 1
11 9 15 1
9 11 12 1
9 11 13 1
9 11 14 1
12 11 13 1
12 11 14 1
13 11 14 1
9 15 16 1
9 15 17 1
16 15 17 1
15 17 18 1
15 17 19 1
18 17 19 1
17 19 20 1
17 19 21 1
17 19 22 1
20 19 21 1
20 19 22 1
21 19 22 1
[ dihedrals ]
; ai aj ak al funct c0 c1 c2 c3 c4 c5
2 1 5 6 9
2 1 5 7 9
3 1 5 6 9
3 1 5 7 9
4 1 5 6 9
4 1 5 7 9
1 5 7 8 9
1 5 7 9 9
6 5 7 8 9
6 5 7 9 9
5 7 9 10 9
5 7 9 11 9
5 7 9 15 9
8 7 9 10 9
8 7 9 11 9
8 7 9 15 9
7 9 11 12 9
7 9 11 13 9
7 9 11 14 9
10 9 11 12 9
10 9 11 13 9
10 9 11 14 9
15 9 11 12 9
15 9 11 13 9
15 9 11 14 9
7 9 15 16 9
7 9 15 17 9
10 9 15 16 9
10 9 15 17 9
11 9 15 16 9
11 9 15 17 9
9 15 17 18 9
9 15 17 19 9
16 15 17 18 9
16 15 17 19 9
15 17 19 20 9
15 17 19 21 9
15 17 19 22 9
18 17 19 20 9
18 17 19 21 9
18 17 19 22 9
[ dihedrals ]
; ai aj ak al funct c0 c1 c2 c3
1 7 5 6 4
5 9 7 8 4
9 17 15 16 4
15 19 17 18 4
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include water topology
#include "amber99sb-ildn.ff/tip3p.itp"
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif
; Include topology for ions
#include "amber99sb-ildn.ff/ions.itp"
[ system ]
; Name
Protein in water
[ molecules ]
; Compound #mols
Protein 1
SOL 876
2. 不同力场可能要求不同的原子名称#
本库AMBER99SB-ILDN的NME甲基碳名称为 CH3，AMBER14SB要求对应名称为 C。资料包中提供两个已经匹配的输入文件；这里仅调整名称，不移动初始几何。
gmx pdb2gmx -f alanine-dipeptide-amber14.pdb -o processed.gro \
-p topol.top -ff amber14sb -water opc -ignh
选择AMBER14SB时使用本篇2026.4模块。该输入库有自己的原子命名、参数及水模型配套文件，不应通过修改共享力场库或强行改残基名称来跳过不匹配。
资料包中的对照流程：
bash run.sh cap99-tip3p
bash run.sh cap14-tip3p
bash run.sh cap14-opc
前两组展示两套蛋白参数库配TIP3P；第三组同时改变水模型为OPC。不要把第三组的差异只归因于某一个因素。
3. 用普通数字输入选择端基#
本例OPLS-AA中，Ala-Ala默认NH3⁺/COO⁻端基总电荷为0。要指定NH3⁺/COOH、形成+1状态，使用 -ter，在本版本菜单输入0和2：
printf '0\n2\n' | gmx pdb2gmx -f ala-ala.pdb \
-o processed.gro -p topol.top -ff oplsaa -water tip3p -ignh -ter
原生菜单里对应 NH3+ 与 COOH。后续必须中和，完整流程是：
bash run.sh alaala-opls-positive
bash run.sh alaala-opls-positive-salt
单残基丙氨酸则需要同时适合两端的模板。本例用OPLS-AA明确标注“仅用于单残基”的两性离子选项，输入1和1：
printf '1\n1\n' | gmx pdb2gmx -f alanine-single.pdb \
-o processed.gro -p topol.top -ff oplsaa -water tip3p -ignh -ter
不要把“仅用于单残基”的选项套给两个残基的肽。不同版本菜单编号可能变化，应核对文字标签；生成后的总电荷必须与所选分子状态一致。这里指定质子化状态，不用这段短MD推断pKa或真实pH下的占比。
4. 运行之后怎样提取phi、psi#
本包为每个封端模型按实际原子顺序生成二面角索引。以AMBER99SB-ILDN/TIP3P分支为例：
[ phi ]
5 7 9 15
[ psi ]
7 9 15 17
进入该分支目录，用原生工具提取：
printf 'phi\n' | gmx angle -f md.trr -n angles.ndx \
-type dihedral -ov phi.xvg
printf 'psi\n' | gmx angle -f md.trr -n angles.ndx \
-type dihedral -ov psi.xvg
phi对应前一残基C—N—CA—C，psi对应N—CA—C—后一残基N。原子序号要跟实际力场生成的顺序一致，不能随便照抄别的体系或另一套参数生成的索引。
bash analyze.sh
数据在 reference/analysis/条件名-dihedral.csv，脚本同时从原始TRR独立重建二面角并与 gmx angle 对照。图中的20ps用于演示提取流程，不表示构象分布或自由能已经收敛；计算自由能请继续参考伞形/WHAM和AWH专题。
下载、提交与结果位置#
下载全部输入、真实轨迹、拓扑、能量表和五张图的绘图脚本。解压进入资料包目录，先查看现成结果：
bash analyze.sh
图片在 reference/figures，总表在 reference/analysis/conditions.csv，各个实际条件的完整文件在 reference/条件名。
从头运行一个水盒示例：
sbatch run.slurm
squeue -u "$USER"
run.slurm 默认选择 water-tip3p。需要其他条件时，修改其中的 CASE_NAME 一行；所有条件名见 CASES.json。新结果写入 work/条件名。完整18条件对照可在新解压目录执行 bash submit.sh；只生成脚本先看、不提交，用 bash submit.sh --dry-run。已有工作目录不会被覆盖。
模板默认核时分区 gpu-8p100-hs，一张 P100、两个 CPU 核；套餐用户按自己可用的分区修改 #SBATCH --partition，例如 P100 套餐分区 gpu-8p100-tc。正文中的单条 bash run.sh ... 仅在已分配的任务资源里使用，不在登录/管理节点运行 MD。
统一环境：gromacs/2026.4-p100-threadmpi，GROMACS2026.4、GNU13.3、CUDA12.8、FFTW3.3.10、thread-MPI/OpenMP。建模、预处理和最小化使用 CPU；动力学阶段短程非键和 PME 使用 GPU，更新/约束使用 CPU。提取使用原生 GROMACS，数值与结构处理使用 NumPy、MDAnalysis2.9.0，图片由 Matplotlib 生成。
本组每个条件都从建模开始，完成最小化、5ps NVT、10ps NPT和20ps生产NVT。它们是短时操作实例；固定0.9nm截断、DispCorr=no和短平衡不用于给水模型、力场精度或性能排名。正式研究按所选参数集文献设置相容的截断/长程处理，并充分平衡与采样。