如题，记录一次分子动力学模拟。体系为 2 mol/L 三氟甲烷磺酸锌 (Zn(OTF)2)，其中加入少量碘化胆碱(ChI)。
使用软件为 gromacs 2024.2 GPU版本，在一块笔记本 GTX 1660Ti 上完成全部计算。
本文全程使用 Amber 力场。读者若想用别的力场也可以，OPLS-AA 和 CHARMM 等都是论文中很常用的力场，模拟有机分子和简单离子的水溶液都比较靠谱。
目录（点击直达，右下角提供快速回顶按钮）
- 目录（点击直达，右下角提供快速回顶按钮）
- 1. 准备工作：力场参数和初始结构 (AmberTools + ACPYPE)
- 2. 计算分子数量和盒子尺寸
- 3. 使用Packmol构建初始盒子
- 4. 创建GROMACS拓扑文件
topol.top - 5. 转换坐标并定义盒子 (gmx editconf)
- 6. 能量最小化
- 7. NVT (恒温恒容) 平衡
- 8. NPT (恒温恒压) 平衡
- 9. 生产阶段模拟
- 10. 数据分析：RDF和配位数
- 11. 总结与其他
1. 准备工作：力场参数和初始结构 (AmberTools + ACPYPE)
这一步至关重要，靠谱的模拟基于正确的力场参数。对于标准离子和水，我们可以从现有力场中提取。对于非标准分子（如胆碱阳离子 CHL+ 和三氟甲磺酸根阴离子 OTF-），我们需要借助外部工具生成参数；或者有文献提供参数的话也可以拿来用。下面我展示如何通过工具获得靠谱的力场参数。
1.1 获取胆碱(CHL)和三氟甲磺酸根(OTF)的力场参数 (AmberTools + ACPYPE)
CHL 和 OTF 不是蛋白质或核酸的标准残基，因此需要专门为其生成 General Amber Force Field（GAFF/GAFF2）参数。这里将详细展示如何通过 AmberTools 套件（antechamber, parmchk2, tleap）为目标分子（以OTF-为例）生成Amber格式的力场参数，然后使用 ACPYPE 将其转换为 GROMACS 格式。CHL+ 的参数化过程与此类似。
注意，既然接触计算模拟，我就默认读者有 Linux 前置技能，掌握简单工具的安装方法。AmberTools 官网 为最新版提供了两种安装方法，拉取conda环境或直接编译二进制都可以，这里建议用conda。ACPYPE 可以通过 pip install acpype 安装。
具体步骤 (以OTF-为例)：
1.1.1 准备初始3D结构
你需要 CHL+ 和 OTF- 的单分子3D结构文件，通常是 .pdb 格式。这些可以通过化学建模软件（如Avogadro、GaussView、Materials Studio）构建并进行初步的几何优化得到（用 DFT 计算优化后的更好），或者从化学数据库下载。确保原子命名规范。
- 例如，我们有一个名为
OTF_initial.pdb的文件作为OTF-的初始结构。
1.1.2 使用antechamber生成MOL2文件和初步原子类型
antechamber是AmberTools中的一个关键程序，用于识别原子类型、分配键类型，并计算原子点电荷。
# 为OTF⁻ (三氟甲磺酸根) 生成MOL2文件
antechamber -i OTF_initial.pdb -fi pdb -o OTF.mol2 -fo mol2 -c bcc -nc -1 -rn OTF -at gaff2
-i OTF_initial.pdb: 输入的PDB文件名。-fi pdb: 输入文件格式为PDB。-o OTF.mol2: 输出MOL2文件名。MOL2文件包含了原子坐标、连接性以及antechamber分配的GAFF2原子类型。-fo mol2: 输出文件格式为MOL2。-c bcc: 使用AM1-BCC方法计算原子点电荷。这是一种快速且对于凝聚相模拟效果较好的电荷计算方法。-nc -1: 指定分子的净电荷为-1。对于CHL+，这里应为-nc 1。-rn OTF: 指定残基名为OTF。-at gaff2: 指定使用GAFF2原子类型。如果想用GAFF，则为-at gaff。
此命令会生成 OTF.mol2 文件。务必检查该文件中各原子的电荷之和是否等于指定的净电荷（这里是 -1）。
1.1.3 使用parmchk2检查并生成缺失的力场参数
GAFF2力场虽然通用，但对于某些特定的键、角、二面角组合可能没有预设参数。parmchk2程序会检查MOL2文件中的原子类型，并与力场库比较，为缺失的参数提供合理的估计值。
parmchk2 -i OTF.mol2 -f mol2 -o OTF.frcmod -s gaff2
-i OTF.mol2: 输入的MOL2文件名 (由antechamber生成)。-f mol2: 输入文件格式为MOL2。-o OTF.frcmod: 输出frcmod (force field modification) 文件名。此文件包含了antechamber无法从标准GAFF2库中找到的所有参数。-s gaff2: 指定参考的力场为GAFF2 (与antechamber中使用的对应)。
打开 OTF.frcmod 文件查看是否有 “ATTN, need revision” 这样的警告，这表示某些参数的估计可能不太可靠，需要进一步检查。如果里面只有一些[ xxx ]标签，下面不带参数，说明所有力场参数全部找到，没有问题。
本文涉及的分子在这一步不会出现问题，都是些很常见的元素。
1.1.4 (可选但推荐) 使用tleap生成Amber的prmtop和inpcrd文件
tleap是Amber中用于整合参数并构建拓扑 (.prmtop) 和坐标 (.inpcrd) 文件的工具。使用tleap生成的prmtop文件可以更好地被ACPYPE处理，因为它已经包含了所有力场参数。
首先，创建一个 tleap 输入脚本，例如 tleap_OTF.in:
source leaprc.gaff2 # 加载GAFF2力场 (或 leaprc.gaff)
loadamberparams OTF.frcmod # 加载parmchk2生成的缺失参数
OTF = loadmol2 OTF.mol2 # 加载分子的MOL2文件
check OTF # 检查分子结构，确保没有错误
saveamberparm OTF OTF.prmtop OTF.inpcrd # 保存Amber拓扑和坐标文件
quit
然后运行 tleap:
tleap -s -f tleap_OTF.in > tleap_OTF.out
- 这会生成
OTF.prmtop(Amber拓扑文件) 和OTF.inpcrd(Amber坐标文件)。 - 检查
tleap_OTF.out文件是否有错误或警告信息。
1.1.5 使用acpype转换为GROMACS格式
ACPYPE (AnteChamber PYthon Parser interface) 可以将Amber的prmtop和inpcrd文件转换为GROMACS兼容的.itp和.gro文件。
acpype -p OTF.prmtop -x OTF.inpcrd -b OTF -o gmx -d
-p OTF.prmtop: 输入Amber拓扑文件。-x OTF.inpcrd: 输入Amber坐标文件。-b OTF: 指定输出文件的前缀 (base name)，例如生成OTF_GMX.itp。-o gmx: 指定输出GROMACS格式。-d(可选): 运行在调试模式，输出更多信息。
此命令会在名为 OTF.acpype 的子目录中生成 OTF_GMX.itp 和 OTF_GMX.gro。
注意，据我测试，最新的2025版会生成 xxx_GMX.top 和 xxx_GMX.gro，你需要根据自己的实际输出来调整。我们需要的 .itp 其实就是所生成的 .top 文件，把.top重命名为.itp即可。核心是获取包含 [ moleculetype ] 定义的GROMACS拓扑片段。
对CHL+也重复以上1.1.1到1.1.5的步骤。
最终产物：
CHL.itp,CHL.groOTF.itp,OTF.gro
恭喜，你已经完全玩明白了：如何获得任一有机分子的力场参数。相信任何小分子溶液的模拟都难不倒你。
1.2 获取Zn2+, I- 和 水分子(SOL)的力场参数 (Amber ff99SB-ILDN)
对于标准离子和水模型，我们可以从 GROMACS 内置的力场或已发表的参数中提取。这里我们使用与 Amber ff99SB-ILDN 兼容的参数。
- Zn2+ (ZN) 和 I- (Iodide, I):
- 这些是单原子离子。它们的
.itp文件非常简单，主要包含[ moleculetype ]和[ atoms ]部分。非键参数（Lennard-Jones σ 和 ε）将放在atomtypes.itp中。 ZN.itp:[ moleculetype ] ; name nrexcl ZN 1 [ atoms ] ; nr type resnr residue atom cgnr charge mass 1 ZN 1 ZN ZN 1 2.0 65.38I.itp:[ moleculetype ] ; name nrexcl I 1 [ atoms ] ; nr type resnr residue atom cgnr charge mass 1 I 1 I I 1 -1.0 126.90
- 这些是单原子离子。它们的
- 水 (SOL):
- 我们使用 tip3p 水模型。其
.itp文件定义了原子、电荷和[ settles ]指令（用于约束水分子几何）。 SOL.itp:[ moleculetype ] ; molname nrexcl SOL 2 [ atoms ] ; nr type resnr residu atom cgnr charge mass 1 OW 1 SOL OW 1 -0.83400 15.99940 2 HW 1 SOL HW1 1 0.41700 1.00800 3 HW 1 SOL HW2 1 0.41700 1.00800 [ settles ] ; OW funct dOH dHH 1 1 0.09572 0.15139 [ exclusions ] 1 2 3 2 1 3 3 1 2
- 我们使用 tip3p 水模型。其
1.3 创建单分子PDB/GRO文件
除了 .itp 文件，我们还需要每个分子类型的单分子坐标文件，作为 Packmol 构建初始盒子的输入。
CHL.pdb和OTF.pdb：可以直接使用 ACPYPE 步骤 (1.1.5) 中生成的.gro文件转换得到 (例如gmx editconf -f OTF.acpype/OTF_GMX.gro -o OTF.pdb)，当然更快的是使用原始的优化过的PDB（确保与ITP中的原子顺序和命名一致）。ZN.pdb:ATOM 1 ZN ZN A 1 0.000 0.000 0.000 1.00 0.00 ZN TER ENDI.pdb:ATOM 1 I I A 1 0.000 0.000 0.000 1.00 0.00 I TER ENDSOL.pdb:ATOM 1 OW SOL A 1 0.000 0.000 0.000 1.00 0.00 O ATOM 2 HW1 SOL A 1 0.096 0.000 0.000 1.00 0.00 H ATOM 3 HW2 SOL A 1 -0.024 0.093 0.000 1.00 0.00 H TER END(注意：水的实际坐标不重要，只要原子名和残基名正确即可。上面是一个示例。)
1.4 整合原子类型参数到 atomtypes.itp
GROMACS需要一个文件来定义所有原子类型的非键相互作用参数（Lennard-Jones σ 和 ε）。我们将所有这些参数集中到 atomtypes.itp 文件中（个人认为是个好习惯，预防后续乱七八糟的报错）。
- 来自 CHL.itp 和 OTF.itp (通过AmberTools/ACPYPE生成): 这些
.itp文件中定义的原子类型 (例如c3,n4,s6,f等) 其对应的 Lennard-Jones 参数源自 GAFF2。你需要从 ACPYPE 生成的 GROMACS 拓扑文件或相关的 Amber 参数文件中提取这些LJ参数（σ 和 ε）。通常，ACPYPE 生成的.itp文件会在其注释部分包含这些参数，或者你可以从prmtop文件解码得到（例如使用parmed库）。 - 来自 Amber ff99SB-ILDN (或其他兼容力场) 的 ZN, I, OW, HW 参数: 这些参数通常在力场的
ffnonbonded.itp文件中。
atomtypes.itp 文件内容示例：
[ atomtypes ]
;name bond_type mass charge ptype sigma epsilon ; comments
; From CHL.itp (GAFF2-like, AmberTools/ACPYPE)
n4 n4 0.00000 0.00000 A 2.49951e-01 1.62122e+01 ; Example N in Choline (CHECK YOUR GAFF2 VALUES)
c3 c3 0.00000 0.00000 A 3.39771e-01 4.51035e-01 ; Example C in Choline & OTF (CHECK YOUR GAFF2 VALUES)
oh oh 0.00000 0.00000 A 3.24287e-01 3.89112e-01 ; Example O in Choline alcohol (CHECK YOUR GAFF2 VALUES)
hx hx 0.00000 0.00000 A 1.88746e-01 8.70272e-02 ; Example H attached to C (CHECK YOUR GAFF2 VALUES)
h1 h1 0.00000 0.00000 A 2.42200e-01 8.70272e-02 ; Example H attached to C next to N/O (CHECK YOUR GAFF2 VALUES)
ho ho 0.00000 0.00000 A 5.37925e-02 1.96648e-02 ; Example H in alcohol (CHECK YOUR GAFF2 VALUES)
; From OTF.itp (GAFF2-like, AmberTools/ACPYPE)
s6 s6 0.00000 0.00000 A 4.05840e-01 2.56898e-01 ; Example S in OTF (CHECK YOUR GAFF2 VALUES)
o o 0.00000 0.00000 A 3.04812e-01 6.12119e-01 ; Example O in OTF (SO3 part) (CHECK YOUR GAFF2 VALUES)
f f 0.00000 0.00000 A 3.03422e-01 3.48109e-01 ; Example F in OTF (CHECK YOUR GAFF2 VALUES)
; For SOL (LJ parameters from TIP3P Amber used with SPC/E charges/geometry)
OW OW 0.00000 0.00000 A 0.315061 0.636386 ; OWT3 from Amber TIP3P
HW HW 0.00000 0.00000 A 0.000000 0.000000 ; HWT3 from Amber TIP3P (LJ usually zero)
; For ZN ion (from amber99sb-ildn.ff/ffnonbonded.itp or compatible source)
ZN ZN 0.00000 0.00000 A 1.95998e-01 5.23000e-02
; For I ion (from amber99sb-ildn.ff/ffnonbonded.itp or compatible source)
I I 0.00000 0.00000 A 4.18722e-01 1.67360e+00
注意：mass 和 charge 列在这里通常设为0，因为实际的质量和电荷在各个分子的 .itp 文件中的 [ atoms ] 部分定义。ptype A 表示普通原子。sigma (nm) 和 epsilon (kJ/mol) 是Lennard-Jones参数。务必确保 atomtypes.itp 中GAFF2相关的原子类型及其参数与你通过AmberTools流程生成的参数一致。
2. 计算分子数量和盒子尺寸
根据体系浓度和期望的模拟规模确定各种分子的数量。
- Zn(OTF)2 浓度: 2 mol/L
- ChI: 31.3 mg (分子量 Choline+ ~104.17 g/mol, I- 126.90 g/mol, ChI ~231.07 g/mol)
- ChI 摩尔数 = 0.0313 g / 231.07 g/mol ≈ 0.00013545 mol
- 宏观 [ChI] = 0.00013545 mol / 0.01 L = 0.013545 mol/L
我们选择在模拟盒子中固定以下数量：
- 胆碱阳离子 (CHL): 1 个
- 碘离子 (I): 1 个
- 锌离子 (ZN): 50 个
- 三氟甲磺酸根离子 (OTF): 100 个 (因为 Zn(OTF)2, 所以 OTF 是 Zn 的两倍)
电荷平衡检查: CHL(+1) + I(-1) + ZN(50 * +2) + OTF(100 * -1) = 1 - 1 + 100 - 100 = 0. 系统电中性。
估算盒子尺寸和水分子数量： 基于50个Zn2+离子和2 mol/L的目标浓度：
- V_box = N_ZN / ([Zn(OTF)2] * N_A)
- V_box = 50 / (2 mol/L * 6.022e23 mol-1) = 4.151e-23 L = 41.51 nm3
- 立方盒子边长
L= (41.51 nm3)^(1/3) ≈ 3.46 nm (34.6 Å)
估算溶质体积 V_solutes ≈ 11.55 nm3。
可用于水的体积 V_water_space = V_box - V_solutes ≈ 29.96 nm3。
每个水分子的体积约 0.03 nm3。
- 水分子 (SOL) 数量
N_SOL≈ 29.96 / 0.03 ≈ 999 个。
最终分子数量：
- CHL: 1
- I: 1
- ZN: 50
- OTF: 100
- SOL: 999
3. 使用Packmol构建初始盒子
Packmol是一个强大的工具，用于将不同类型的分子按照指定的数量和空间约束填充到模拟盒子中。
packmol.inp 文件：
# Packmol input file for Zn(OTF)2 / ChI aqueous solution
tolerance 2.0
output system_packed.pdb
filetype pdb
# Box dimensions: 3.46 nm = 34.6 Angstroms
structure CHL.pdb
number 1
inside cube 0. 0. 0. 34.6
end structure
structure I.pdb
number 1
inside cube 0. 0. 0. 34.6
end structure
structure ZN.pdb
number 50
inside cube 0. 0. 0. 34.6
end structure
structure OTF.pdb
number 100
inside cube 0. 0. 0. 34.6
end structure
structure SOL.pdb
number 999
inside cube 0. 0. 0. 34.6
end structure
其中，tolerance 指每个插入的分子之间的最小距离，一般设为 2.0 就可以。
运行 Packmol:
packmol < packmol.inp
这将生成 system_packed.pdb 文件。
4. 创建GROMACS拓扑文件 topol.top
topol.top 文件是GROMACS模拟的核心，它定义了系统的力场、包含的分子类型及其数量。
topol.top 文件内容：
; Topology file for Zn(OTF)2 / ChI aqueous solution
; 1. Force field parameters
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 2 yes 0.5 0.8333 ; AMBER/GAFF standard
; Include the non-bonded atom type definitions
#include "atomtypes.itp"
; 2. Include molecule topology (.itp files)
#include "CHL.itp" ; Choline cation
#include "I.itp" ; Iodide anion
#include "ZN.itp" ; Zinc cation
#include "OTF.itp" ; Triflate anion
#include "SOL.itp" ; Water
; 3. System definition
[ system ]
; Name of the system
Zn(OTF)2 and ChI in Water
; 4. Molecules
; List of molecule types and their counts in the system.
[ molecules ]
; Compound nmol
CHL 1
I 1
ZN 50
OTF 100
SOL 999
重要：确保 topol.top、所有 .itp 文件和 system_packed.pdb 在同一工作目录。
至于topol文件中应该包含哪些东西，注释里已经写得很清楚了，用#include把已放进同一目录中，所有该引用的文件全引进来。最后的分子数和前面的步骤对应，别搞错了。
注意：在 topol.top 的 #include "atomtypes.itp" 之前记得添加 [ defaults ] ，用于正确初始化力场设置，否则会报语法错误。
5. 转换坐标并定义盒子 (gmx editconf)
将Packmol生成的PDB文件转换为GROMACS的GRO格式，并确保盒子尺寸正确。
gmx editconf -f system_packed.pdb -o system.gro -c -box 3.46 3.46 3.46
-f system_packed.pdb: 输入PDB。-o system.gro: 输出GRO。-c: 将分子置于盒子中心。-box 3.46 3.46 3.46: 显式定义盒子尺寸 (单位nm)。
检查 system.gro 文件最后一行的盒子向量和第二行的总原子数 (应为 1*21 + 1*1 + 50*1 + 100*8 + 999*3 = 3869 个原子)。
6. 能量最小化
移除初始结构中的不良接触，降低系统能量。
em.mdp 文件：
; MDP file for Energy Minimization
define = -DFLEXIBLE
integrator = steep
emtol = 1000.0
emstep = 0.01
nsteps = 50000
nstlist = 10
cutoff-scheme = Verlet
pbc = xyz
coulombtype = PME
rcoulomb = 1.0
vdwtype = Cut-off
rvdw = 1.0
DispCorr = EnerPres
constraints = h-bonds
constraint_algorithm = LINCS
运行Grompp和MDrun：
gmx grompp -f em.mdp -c system.gro -p topol.top -o em.tpr -maxwarn 1
gmx mdrun -v -deffnm em
检查：em.log是否显示 “Steepest Descents converged to Fmax < 1000”。我这里在576步收敛，Fmax ≈ 957 kJ/mol/nm，最终势能 ≈ -8.9e4 kJ/mol。
7. NVT (恒温恒容) 平衡
在恒定体积和目标温度下弛豫系统。
nvt.mdp 文件：
; MDP file for NVT equilibration
title = ZnOTF2_ChI NVT equilibration
define = -DFLEXIBLE
integrator = md
nsteps = 50000 ; 100 ps
dt = 0.002
nstlog = 5000
nstenergy = 5000
nstxout-compressed = 5000
continuation = no
constraint_algorithm = lincs
constraints = h-bonds
cutoff-scheme = Verlet
nstlist = 20
pbc = xyz
coulombtype = PME
rcoulomb = 1.0
vdwtype = Cut-off
rvdw = 1.0
DispCorr = EnerPres
tcoupl = V-rescale
tc-grps = System
tau_t = 0.1
ref_t = 298.15 ; Target temperature (K)
pcoupl = no
gen_vel = yes
gen_temp = 298.15
gen_seed = -1
运行Grompp和MDrun：
gmx grompp -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr -r em.gro -maxwarn 1
gmx mdrun -v -deffnm nvt
检查：使用 gmx energy -f nvt.edr -o temperature_nvt.xvg (选择 “Temperature”)。
我这里显示平均温度为 298.065 K，接近目标值298.15 K，波动合理。
8. NPT (恒温恒压) 平衡
在恒定压力和目标温度下弛豫系统，允许盒子体积变化以达到正确的密度。
npt.mdp 文件：
; MDP file for NPT equilibration
title = ZnOTF2_ChI NPT equilibration
define = -DFLEXIBLE
integrator = md
nsteps = 500000 ; 1 ns
dt = 0.002
nstlog = 5000
nstenergy = 5000
nstxout-compressed = 5000
continuation = yes
constraint_algorithm = lincs
constraints = h-bonds
cutoff-scheme = Verlet
nstlist = 20
pbc = xyz
coulombtype = PME
rcoulomb = 1.0
vdwtype = Cut-off
rvdw = 1.0
DispCorr = EnerPres
tcoupl = V-rescale
tc-grps = System
tau_t = 0.1
ref_t = 298.15
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 2.0
ref_p = 1.0 ; Target pressure (bar)
compressibility = 4.5e-5
refcoord_scaling = com
运行Grompp和MDrun：
gmx grompp -f npt.mdp -c nvt.gro -p topol.top -o npt.tpr -t nvt.cpt -r nvt.gro -maxwarn 1
gmx mdrun -v -deffnm npt
检查：使用 gmx energy 分析 npt.edr 文件中的 Temperature, Pressure, Density (或 Volume)。确保它们都已达到稳定。
9. 生产阶段模拟
系统充分平衡后（即温度、压力、密度等关键物理量在目标值附近稳定波动），就可以进行长时间的生产阶段分子动力学模拟了。这个阶段的主要目的是采集足够的数据样本，用于后续的各种物理性质分析和现象观察。
md.mdp 文件：
; MDP file for Production MD
title = ZnOTF2_ChI Production MD ; 描述性标题
define = -DFLEXIBLE ; 与平衡阶段一致，通常用于水模型
; Run parameters
integrator = md ; leap-frog integrator, 标准的MD积分器
nsteps = 5000000 ; 总模拟步数。这里是 5,000,000 步
; 5,000,000 * 0.002 ps/step = 10,000 ps = 10 ns
; **这是需要根据研究目的和计算资源仔细设定的关键参数。**
; 对于离子溶液的结构和动力学，通常需要几十到几百纳秒。
dt = 0.002 ; 时间步长 (ps)。2 fs是常见选择，前提是约束了含氢键。
; Output control - 控制各种输出文件的频率
nstxout = 0 ; 输出全精度坐标到 .trr 文件的频率。设为0表示不输出。
; 通常在生产阶段不保存 .trr，因为 .xtc 已足够且小得多。
nstvout = 0 ; 输出速度到 .trr 文件的频率。设为0表示不输出。
nstfout = 0 ; 输出力到 .trr 文件的频率。设为0表示不输出。
nstlog = 10000 ; 每10000步 (20 ps) 更新一次 .log 日志文件。
; 对于长模拟，可以适当增大此值以减小日志文件大小。
nstenergy = 10000 ; 每10000步 (20 ps) 保存一次能量项到 .edr 文件。
; 用于后续热力学性质分析。
nstxout-compressed = 10000 ; 每10000步 (20 ps) 保存一次压缩坐标到 .xtc 文件。
; **这是用于后续轨迹分析的主要文件。**
; 频率决定了分析的时间分辨率。
; Bond parameters
continuation = yes ; **关键！** 表示从上一个模拟阶段 (NPT) 继续。
; 会读取上一步的检查点文件 (.cpt) 以获取坐标、速度、盒子状态等。
constraint_algorithm = lincs ; LINear Constraint Solver，用于约束键长。
constraints = h-bonds ; 约束所有包含氢原子的键。这允许使用更大的时间步长(如2fs)。
lincs_iter = 1 ; LINCS算法的迭代次数，影响约束精度。
lincs_order = 4 ; LINCS算法的阶数，也影响精度。
; Neighbor searching - 邻近列表参数，与平衡阶段一致
cutoff-scheme = Verlet
nstlist = 20 ; 每20步 (40 fs) 更新一次邻近列表。
pbc = xyz ; 在所有三个方向上使用周期性边界条件。
; Electrostatics - 静电相互作用参数，与平衡阶段一致
coulombtype = PME ; Particle Mesh Ewald，用于处理长程静电作用。
rcoulomb = 1.0 ; 短程静电截断半径 (nm)。
; Van der Waals - 范德华相互作用参数，与平衡阶段一致
vdwtype = Cut-off ; 使用截断处理范德华相互作用。
rvdw = 1.0 ; 范德华截断半径 (nm)。
DispCorr = EnerPres ; 对能量和压力进行长程色散校正。
; Temperature coupling - 温度耦合参数，与平衡阶段一致
tcoupl = V-rescale ; V-rescale 恒温器，能较好地产生正则系综。
tc-grps = System ; 将整个系统耦合到同一个热浴。
tau_t = 0.1 ; 温度耦合的时间常数 (ps)。
ref_t = 298.15 ; 目标参考温度 (K)。
; Pressure coupling - 压力耦合参数，与平衡阶段一致
pcoupl = Parrinello-Rahman ; Parrinello-Rahman 恒压器，允许盒子形状和体积变化。
pcoupltype = isotropic ; 各向同性压力耦合，盒子在所有方向上等比例缩放。
tau_p = 2.0 ; 压力耦合的时间常数 (ps)。
ref_p = 1.0 ; 目标参考压力 (bar)。
compressibility = 4.5e-5 ; 体系的压缩系数 (bar⁻¹)，这里使用水的典型值。
refcoord_scaling = com ; 参考坐标缩放方式，对没有位置限制的体系影响不大。
; Velocity generation (NOT needed if continuation = yes)
; gen_vel = no ; 由于 continuation = yes，不需要重新生成速度。```mdp
; MDP file for Production MD
title = ZnOTF2_ChI Production MD ; 描述性标题
define = -DFLEXIBLE ; 与平衡阶段一致，通常用于水模型
; Run parameters
integrator = md ; leap-frog integrator, 标准的MD积分器
nsteps = 5000000 ; 总模拟步数。这里是 5,000,000 步
; 5,000,000 * 0.002 ps/step = 10,000 ps = 10 ns
; **这是需要根据研究目的和计算资源仔细设定的关键参数。**
; 对于离子溶液的结构和动力学，通常需要几十到几百纳秒。
dt = 0.002 ; 时间步长 (ps)。2 fs是常见选择，前提是约束了含氢键。
; Output control - 控制各种输出文件的频率
nstxout = 0 ; 输出全精度坐标到 .trr 文件的频率。设为0表示不输出。
; 通常在生产阶段不保存 .trr，因为 .xtc 已足够且小得多。
nstvout = 0 ; 输出速度到 .trr 文件的频率。设为0表示不输出。
nstfout = 0 ; 输出力到 .trr 文件的频率。设为0表示不输出。
nstlog = 10000 ; 每10000步 (20 ps) 更新一次 .log 日志文件。
; 对于长模拟，可以适当增大此值以减小日志文件大小。
nstenergy = 10000 ; 每10000步 (20 ps) 保存一次能量项到 .edr 文件。
; 用于后续热力学性质分析。
nstxout-compressed = 10000 ; 每10000步 (20 ps) 保存一次压缩坐标到 .xtc 文件。
; **这是用于后续轨迹分析的主要文件。**
; 频率决定了分析的时间分辨率。
; Bond parameters
continuation = yes ; **关键！** 表示从上一个模拟阶段 (NPT) 继续。
; 会读取上一步的检查点文件 (.cpt) 以获取坐标、速度、盒子状态等。
constraint_algorithm = lincs ; LINear Constraint Solver，用于约束键长。
constraints = h-bonds ; 约束所有包含氢原子的键。这允许使用更大的时间步长(如2fs)。
lincs_iter = 1 ; LINCS算法的迭代次数，影响约束精度。
lincs_order = 4 ; LINCS算法的阶数，也影响精度。
; Neighbor searching - 邻近列表参数，与平衡阶段一致
cutoff-scheme = Verlet
nstlist = 20 ; 每20步 (40 fs) 更新一次邻近列表。
pbc = xyz ; 在所有三个方向上使用周期性边界条件。
; Electrostatics - 静电相互作用参数，与平衡阶段一致
coulombtype = PME ; Particle Mesh Ewald，用于处理长程静电作用。
rcoulomb = 1.0 ; 短程静电截断半径 (nm)。
; Van der Waals - 范德华相互作用参数，与平衡阶段一致
vdwtype = Cut-off ; 使用截断处理范德华相互作用。
rvdw = 1.0 ; 范德华截断半径 (nm)。
DispCorr = EnerPres ; 对能量和压力进行长程色散校正。
; Temperature coupling - 温度耦合参数，与平衡阶段一致
tcoupl = V-rescale ; V-rescale 恒温器，能较好地产生正则系综。
tc-grps = System ; 将整个系统耦合到同一个热浴。
tau_t = 0.1 ; 温度耦合的时间常数 (ps)。
ref_t = 298.15 ; 目标参考温度 (K)。
; Pressure coupling - 压力耦合参数，与平衡阶段一致
pcoupl = Parrinello-Rahman ; Parrinello-Rahman 恒压器，允许盒子形状和体积变化。
pcoupltype = isotropic ; 各向同性压力耦合，盒子在所有方向上等比例缩放。
tau_p = 2.0 ; 压力耦合的时间常数 (ps)。
ref_p = 1.0 ; 目标参考压力 (bar)。
compressibility = 4.5e-5 ; 体系的压缩系数 (bar⁻¹)，这里使用水的典型值。
refcoord_scaling = com ; 参考坐标缩放方式，对没有位置限制的体系影响不大。
; Velocity generation (NOT needed if continuation = yes)
; gen_vel = no ; 由于 continuation = yes，不需要重新生成速度。
运行Grompp和MDrun：
gmx grompp -f md.mdp -c npt.gro -p topol.top -o md.tpr -t npt.cpt -r npt.gro -maxwarn 1
#-f md.mdp: 指定MDP参数文件。
#-c npt.gro: 提供NPT平衡结束时的坐标作为起始结构。
#-p topol.top: 指定系统的拓扑文件。
#-o md.tpr: 输出的二进制运行文件。
#-t npt.cpt: **关键！** 读取NPT平衡的检查点文件，包含速度、精确盒子状态等，用于无缝续跑。
#-r npt.gro: (可选，但推荐) 参考结构，用于某些类型的分析或限制。
#-maxwarn 1: 允许一个警告通过而不中断grompp。
gmx mdrun -v -deffnm md
#-v: 输出详细运行信息到屏幕。
#-deffnm md: GROMACS会自动寻找md.tpr作为输入，
# 并将所有输出文件 (md.log, md.xtc, md.edr, md.gro, md.cpt等) 以"md"为前缀命名。
关于输出文件：
md.log: 文本日志，约几MB到几十MB。md.edr: 二进制能量文件，约几MB到几十MB。md.xtc: 压缩轨迹文件，是分析的主要数据源。对于10ns、约4000原子的体系、每20ps保存一帧（500帧），大小约几十MB到一百多MB。
好了，如果到这里一切顺利，那么恭喜你完成了分子动力学模拟，现在是喜闻乐见的提数据作图环节。
10. 数据分析：RDF和配位数
这里分析Zn2+离子周围O原子(来自OTF-)和O原子(来自水)的径向分布函数和配位数。
10.1 创建索引文件 (gmx make_ndx)
gmx make_ndx -f md.tpr -o index.ndx
在交互提示符下创建以下组：
- 选择Zn原子:
a ZN(这步操作会创建一个新组，我这里是组11，根据make_ndx输出确认，或者空敲一下回车) - 命名Zn组:
name 11 ZN_ions(将11替换为实际组号，主要防止弄错) - 选择OTF中的氧原子 (原子名O1, O2, O3，残基名OTF):
r OTF & a O1 O2 O3(假设创建组12) - 命名OTF氧原子组:
name 12 O_OTF - 选择水分子的氧原子 (原子名OW，残基名SOL):
r SOL & a OW(假设创建组13) - 命名水氧原子组:
name 13 O_Water - 保存并退出:
q
确认：make_ndx 会列出预定义的组，如：
0 System : 3869 atoms
...
4 ZN : 50 atoms (基于残基名ZN的预定义组)
...
确保你使用 a ZN (按原子名) 创建的新组被正确识别和命名。
不知道怎么看？操作完，仍在make_ndx中空敲一下回车键就可以。
10.2 计算RDF和配位数 (gmx rdf)
- Zn - O-OTF:
gmx rdf -s md.tpr -f md.xtc -n index.ndx -o rdf_zn_ootf.xvg -cn cn_zn_ootf.xvg -ref ZN_ions -sel O_OTF -tu ns -bin 0.002 -rmax 1.0 # -s md.tpr: 提供拓扑和原子信息，用于正确识别原子和计算距离。 # -f md.xtc: 输入生产阶段的轨迹文件。 # -n index.ndx: 指定包含原子组定义的索引文件。 # -o rdf_zn_ootf.xvg: 输出RDF数据到此文件。 # -cn cn_zn_ootf.xvg: 输出累积配位数 (Coordination Number) 数据到此文件。 # -ref ZN_ions: 指定参考组 (即中心原子组)。这里是你之前在index.ndx中定义的ZN_ions组。 # -sel O_OTF: 指定选择组 (围绕参考组分布的原子组)。这里是你定义的O_OTF组。 # -tu ns: 指定输出文件中的时间单位为纳秒 (如果轨迹是以纳秒为单位记录的，这有助于xvg文件的标签)。 # -bin 0.002: 设置RDF计算中距离柱状图的宽度 (bin size) 为0.002 nm。 # 较小的值可以得到更平滑的曲线，但可能需要更多数据或更长的模拟。 # -rmax 1.0: 设置计算RDF的最大距离为1.0 nm。通常关注第一和第二配位层，此范围足够。 # 默认值通常是盒子长度的一半。 - Zn - O-Water:
gmx rdf -s md.tpr -f md.xtc -n index.ndx -o rdf_zn_owater.xvg -cn cn_zn_owater.xvg -ref ZN_ions -sel O_Water -tu ns -bin 0.002 -rmax 1.0 # 同上，只不过这里研究Zn - O-Water
解读输出：
rdf_*.xvg: g(r) vs r。第一个峰表示第一配位层距离，第一个谷定义其边界。cn_*.xvg: 配位数 vs r。在RDF第一个谷对应的r值处的配位数值即为第一配位层的平均配位数。
11. 总结与其他
总之，本文涵盖了从力场参数准备、体系模拟到数据分析的关键环节，并给出了所有输入文件的详细内容，理论上完成 #1 准备工作 后可以顺利复现，并举一反三地推广到其他任何体系。
我认为按照此流程，至少跑绝大多数液相中的MD模拟都不在话下，并能获得论文中常见的 RDF、配位数等数据，和花里胡哨的 “满盒子的分子” 图 (这个百度能找到教程，用VMD就可以画出你想要的)。
后续可以进行的分析包括：
- 离子的扩散系数计算 (MSD)。
- 特定分子构象分析。
- 氢键网络分析。
- 更长时间的模拟以获得更好的统计精度或观察慢过程。
对于Gromacs，我也从一无所知中摸索过来，深知教程/信息的零散过时、完整案例的缺失难寻，和海量报错带来的烦躁崩溃。希望这篇博文能对你有哪怕一丝些微的帮助！