GROMACS进行蛋白小分子复合物分子动力学模拟流程及遇到的问题
当体系非中性时会触发警告信息，在代码末尾添加 -maxwarn 1 参数即可实现预期效果。例如基于我的redocking实验结果，请详细阐述这一现象的具体表现及原因分析。 针对出现的问题进行了相应的解释与推测：如果有不准确请指正；其中黄体显示的部分属于可选方案。本文主要阐述了MD正常运行的方法论框架，在具体实施过程中可以根据实际情况灵活调整研究方案。
目录
一、复合物的预处理：
二、蛋白的准备
三、【重要】准备配体：
四、蛋白和配体的组装：
五、定义盒子
六、添加溶剂
七、添加抗衡离子
八、能量最小化
九、位置限制性NVT
十、位置限制性NPT
十一、正式MD
一、复合物的预处理：
RSCB平台（RCSB PDB: Homepage）提供了共晶2FQ9的数据集，在使用Schrodinger软件进行分子模拟之前进行了初步处理（特别注意，在随后的力场计算中可能会因为蛋白质残基缺失而导致错误信息出现，请确保补充完整的蛋白质残基）。随后进行了优化计算（minimize），生成网格（grid），并在红ocking过程中筛选出与靶标RMSD值小于2.0的构象后进行合并（merge），最终导出结果文件（format为2FQ9.pdb）。建议分别将蛋白结构和配体结构单独保存为2FQ9.pdb和ligS3.mol2格式文件（这样可以避免后续操作中的诸多问题）。尽管这种方法主要针对的是共价抑制剂的研究（目的是探讨非共价作用以及疏水亲核效应），但通过非共价对接的方式也能较好地模拟实际实验中的构象特征。
二、蛋白的准备
grep CRJ 2FQ9.pdb > ligS3.pdb//分离辅因子以生成独立的坐标数据文件，在pdb格式中使用的重要参数是crj。
创建并保存一份2FQ9.pdb文件的副本；随后手动移除该文件中的配体行以及连接信息部分（其中连接信息不重要）
gmx pdb2gmx -f 2FQ9.pdb -o 2FQ9_pre.gro -ignh//将只含有蛋白的pdb文件转换为gro文件
-ignh: Omitting all hydrogen atoms in the input PDB file is specified by this command.
This is because the PDB file was generated by NMR and thus contains hydrogen atoms.
It ensures that all hydrogen atoms follow the GROMOS naming conventions.
To prevent any inconsistencies in the hydrogen atom names.
-ff: 设置使用的力场，请根据自身的体系进行选择.若无需提示即可自行设定，则蛋白复合物建议选用Amber99SB-ILDN 力场（亦有更为高级的 Amber14SB 和 Amber19SB 可供参考）。
-f: 指定需要进行处理的蛋白质结构文件.
-o: 指定一个新生成的GROMACS文件名(也可以是其它的文件类型, 如pdb).
-p: 指定新生成的拓扑文件名（包含所有原子及其相互作用参数），若无，默认输出topol.top。
-water: 指定水模型. 不输入-water会提示选择模型，推荐使用TIP3P水模型.
生成了一组输出文件：包含蛋白的结构信息2FQ9_pre.gro、蛋白的主拓扑信息topol.top以及蛋白的位置限制信息posre.itp。其中 gro 用于存储原子坐标信息、top 用于描述分子骨架连接关系、itp则作为主拓扑信息可用于构建。
三、【重要】准备配体：
./RESP2.sh ligS3.pdb将最新版Multiwfn位于examples\RESP目录内的RESP2.sh文件进行拷贝操作至ligS3.pdb所在的文件夹中（Multiwfn对多种输入格式均提供支持），并借助脚本程序计算得到RESP值为0.5的结果，并生成包含resp信息的chg格式文件
antechamber -i ligS3.pdb -fi pdb -o ligS3_resp.mol2 -fo mol2 -pf y由于Amborigens无法解析chg文件而导致无法在antechamber中设置相应的静电势参数建议将配体的PDB文件转换为mol2格式以供处理如果直接从Schrodinger的输出导出mol2格式则无需进行此转换步骤
手动复制.chg中的电荷信息到.mol2文件中（坐标后面一列）。
在处理含有温度因子ANISOU标记的PDB文件时，在线程分析阶段需要完成特定操作。首先通过grep命令筛选关键词生成仅包含配体结构的新PDB文件。具体步骤是首先使用grep提取配体名称并生成带有ANISOU标记的新文件，在此过程中每个原子占据两个状态（各占半数）。接着通过再次应用grep筛选出不带ANISOU标记的部分，并记录这些状态下所有原子均为偶数行的情况。最后使用sed命令过滤出奇数行为的主要构象特征，并将结果保存至lig.pdb中作为最终输出
（可选）
//也可通过Gaussian直接计算RESP，输出.gesp文件（但更推荐上面方法）。
antechamber -i ligS3.gesp -fi gesp -o ligS3.mol2 -fo mol2 -pf y -c resp//利用Ambertools的amtechamber将.resp文件整合进gesp中，生成.mol2文件。
parmchk2 -i ligS3_resp.mol2 -f mol2 -o ligS3.frcmod//parmchk2 是 parmchk 的增强版, 能够识别并补充分子构型中缺少的 GAFF 参数, 并生成相应的补充参数文件 LigS3.frcmod.
tleap //激活tleap命令
source leaprc.gaff //载入力场
loadamberparams ligS3.frcmod //载入配体参数化文件
lig = loadmol2 ligS3_resp.mol2 //载入配体
check lig //检查配体
saveamberparm lig ligS3.prmtop ligS3.inpcrd //保存amber输入文件,得到了prmtop文件和inpcrd文件（可选）
这部分可以通过一个编译好的leap.in文件执行，内容如下：
source leaprc.gaff
loadamberparams ligS3.frcmod
lig=loadmol2 ligS3_resp.mol2
check lig
saveamberparm lig ligS3.prmtop ligS3.inpcrd
quit tleap -f leap.in //运行leap.in acpype -p lig.prmtop -x lig.inpcrd借助Ambertools的acpype工具对文件进行格式转换（从Ambortage到GROMACS），默认输出文件名为.prmtop及lig.inpcrd中的对应单位名称（例如：CRJ_GMX.gro及CRJ_GMX.top）
创建ligS3.itp文件，并从CRJ_GMX.top中提取[ moleculetype ]-[ system ]部分（不包含系统信息），粘贴到ligS3.itp文件中作为配体的拓扑文件。
gmx genrestr -f CRJ_GMX.gro -o posre_ligS3.itp -fc 1000 1000 1000//限制配体, 利用ligS3.gro文件生成限制posre_ligS3.itp文件。
-f CRJ_GMX.gro: 输入配体的结构文件.
-o posre_ligS3.itp: 输出配体施加位置限制的itp文件.
-fc 1000 1000 1000: 位置限制的力常数(单位: kJ/mol-nm2).
四、蛋白和配体的组装：
gro文件：获取一个蛋白的gro文件（2FQ9.gro），命名为2FQ9_complex.gro后打开查看；然后将ligS3.gro的坐标文件复制到该蛋白坐标系的下方位置；由于增加了分子中原子数量后导致总原子数发生变化，请确保第二行中记录的总原子数数据与新增部分相加计算以获得准确值；
//gro文件格式要求极为严格，在复制过程中需与原有列保持一致。ligand的PDB文件末尾无需复制，请仅关注并复制坐标相关的行即可
将力场的文件夹（Amber99SB-ILDN.ff）移动至与现有gro和topol.top文件位于同一目录下的新文件夹中，并将其命名为Amber99SB-ILDN_modify.
//避免污染源文件，因为定义配体原子需要修改力场文件。
打开ffnonbonded.itp文件所在的文件夹，并将其ligS3.itp文件中的**[ atomtype ]**部分内容附加在末尾处。
//定义配体的原子类型，该文件在力场中include。
3）修改topol.top文件：
用于topol.top文件中对所有力场路径进行替换工作，并对其中每个路径进行详细调整以适应新的计算需求
将ligS3.itp文件放入开头位置的力场引入部分，并在其下方使用分号进行标注以表明注释的位置；同时，在此位置上标记出配体分子对应的拓扑结构信息。
; Include forcefield parameters
#include "Amber99SB-ILDN_modify.ff/forcefield.itp"
; Include ligand topology
#include "ligS3.itp"3.特别提示以下几点：第一，在文件末尾配体拓扑文件和力场/水文件之间导入生成的配体限制itp文件posre_ligS3.itp；第二，请确保所导入的定义词"POSRES_LIG"与nvt.mdp及npt.mdp中的define内容"-DPOSRES_LIG"保持一致
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Ligand position restraints
#ifdef POSRES_LIG
#include "posre_ligS3.itp"
#endif
; Include water topology
#include "Amber99SB-ILDN_modify.ff/tip3p.itp"4.优化相关分子部分, 因为新增了配体分子, 应补充配体名称(需参考acpype生成文件中的名称)及其数量.
[ molecules ]
; Compound #mols
Protein_chain_A 1
CRJ 1在此处建议详细说明如何备份整个文件夹
五**、定义盒子**
gmx editconf -f 2FQ9_complex.gro -o 2FQ9_box.gro -bt dodecahedron -d 1.0//使用editconf模块定义盒子大小。
六**、添加溶剂**
gmx solvate -cp 2FQ9_box.gro -cs spc216.gro -p topol.top -o 2FQ9_solv.gro//采用溶质填充整个区域进行处理后，在运行结束后top文件中新增SOL项
七**、添加抗衡离子**
gmx grompp -f minim.mdp -c 2FQ9_solv.gro -p topol.top -o 2FQ9_ions.tpr//如果体系不为中性会有warning，可在后面加 -maxwarn 1 。
gmx genion -s 2FQ9_ions.tpr -o 2FQ9_ions.gro -p topol.top -neutral//添加离子,默认添加NA和CL，-neutral平衡至中性。
在终端界面中设置第15组（SOL），然后以新增的离子替代原本配置的SOL水分子，并按回车键提交命令。操作后可查看topol.top文件以确认新增内容。
（可选）
gmx genion -s 2FQ9_ions.tpr -o 2FQ9_ions.gro -p topol.top -pname NA -nname CL -nn 6//手动添加离子种类及数量。
八**、能量最小化**
gmx grompp -f minim.mdp -c 2FQ9_ions.gro -p topol.top -o em.tpr//为了提高能量最小化步骤的效果数量可以适当增加收敛性约束的程度从而尽量确保结构达到最佳稳定性水平这样的优化调整有助于减少后续可能出现的问题
gmx mdrun -v -deffnm em//报错可尝试加 -nt 1（具体情况具体分析）。
gmx energy -f em.edr -o potential.xvg//检查势能，可用xmgrace查看。
xmgrace potential.xvg九**、位置限制性NVT**
（可选）
gmx make_ndx -f em.gro -o index.ndx
> 1 | 13
> q在nvt.mdp文件中，当tc_grps属性设为'system'时，则无需执行这一步。原因在于直接将配体与其自身的热浴场耦合是不可行的，这会以防止体系发生过度反应或崩溃的前提下仍建议采取此措施。
//此操作涉及将蛋白（1）与配体（13）归为同一组合并处理它们之间的相互作用关系。由于此次分组数目减少至2组，请对各组合中的相关参数进行详细设置；具体而言，在每一步骤中都需要执行相应的参数优化操作以确保模拟结果的一致性与准确性。
生成过程中的每一份molecular dynamics (MDP)文件都需要按照上述设定进行相应的参数配置；在软件界面中选择"参数设置"选项卡后依次完成相关参数的操作即可。
tc-grps = Protein_CRJ Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 310 310 ; reference temperature, one for each group, in K gmx grompp -f nvt.mdp -c em.gro -p topol.top -r em.gro -o nvt.tpr// -maxwarn -1：选择性地忽略 warnings 但仍需予以关注，在后续任务出现问题时，请从错误提示中寻找原因并采取相应措施。处理 nvt.mdp 文件时，请注意调整其时间步长、时间步数等相关参数（如果出现错误）。请务必注意文件名称的位置限制，并确保 define 的命令与 top 文件中的 include 命令之间存在对应关系。对于这一步骤中的 warning 建议不要轻率地忽略，在遇到问题时应优先解决以避免影响最终结果。
// -r（输入约束由所使用的软件版本号决定，在旧版中并非强制性要求；通常与-c传递gro文件相关）
（可选择性地，在NPT和正式MD中若不使用index分组，则遵循无-n index.ndx的对应ndx命令，下同）
gmx grompp -f nvt.mdp -c em.gro -p topol.top -r em.gro -n index.ndx -o nvt.tpr gmx mdrun -v -deffnm nvt十**、位置限制性NPT**
gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -r nvt.gro -o npt.tpr
gmx mdrun -v -deffnm npt十一、正式MD
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr
gmx mdrun -nt 20 -v -deffnm md本任务是在单个GPU并利用OpenMP技术下运行20个线程。由于仅使用了一个GPU的情况，在这种配置下进行多线程处理时这使得线程数量对性能的影响微乎其微。
mdrun的一些命令：
-v可是显示进度信息，并提供预期结束时间 -deffnm所有必需的输出字段默认名称设置为name
-x输出xtc文件
-cpi可以生成备份的参数(.cpt)文件，可以续跑
（可选）
gmx mdrun -pin on -ntmpi 1 -ntomp 12 -pme gpu XXX.tpr//若采用单块高质量GPU，可用GPU计算PME。