GROMACS 蛋白质-配体复合物模拟全流程：从对接结果到轨迹分析

# Ligand prepared with Meeko before docking: restore bond orders and all H from the PDBQT REMARKs
mk_export.py vina_out.pdbqt -s lig_docked.sdf      # older releases use -o; check mk_export.py --help
# Keep only the pose you want (first record = top-ranked pose), then write mol2 for CGenFF / ACPYPE
obabel lig_docked.sdf -O lig.mol2 -l 1
# obabel names the residue UNL1; rename it so ACPYPE writes residue LIG (make_ndx, gmx_MMPBSA and the analysis commands select "LIG")
sed -i.bak 's/UNL1/LIG1/' lig.mol2
# Check the net charge before parameterization
python3 -c "from rdkit import Chem; m=Chem.MolFromMolFile('lig_docked.sdf',removeHs=False); print(Chem.GetFormalCharge(m), m.GetNumAtoms())"

# Ligand NOT prepared with Meeko: rebuild bond orders from a SMILES template, keep the docked coordinates
from rdkit import Chem
from rdkit.Chem import AllChem

pose = Chem.MolFromPDBFile("lig_pose.pdb", removeHs=True)          # heavy atoms of the docked pose
template = Chem.MolFromSmiles("CCCc1ccccc1O")                        # same protonation state you want in MD
pose = AllChem.AssignBondOrdersFromTemplate(template, pose)
pose_h = Chem.AddHs(pose, addCoords=True)                           # H placed on the docked heavy atoms
Chem.MolToMolFile(pose_h, "lig_docked.sdf")
print("formal charge:", Chem.GetFormalCharge(pose_h))
# Then, as above: obabel lig_docked.sdf -O lig.mol2 && sed -i.bak 's/UNL1/LIG1/' lig.mol2

# Protein only: remove ligand, crystal waters and additives first
grep -v -e HETATM -e CONECT complex_docked.pdb > protein.pdb

# CHARMM36m route (force field directory unpacked in the working directory)
# -ter: choose NH3+ and COO-. Without it an N-terminal MET gets MET1 from ethers.n.tdb and pdb2gmx stops
gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top -ff charmm36-feb2026_cgenff-5.0 -water tip3p -ignh -his -ter

# AMBER route (amber14sb.ff ships with GROMACS 2026)
gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top -ff amber14sb -water tip3p -ignh -his

# GAFF2 + AM1-BCC; always pass the net charge with -n
acpype -i lig.mol2 -b LIG -c bcc -a gaff2 -n 0
cd LIG.acpype
# LIG_GMX.itp starts with [ atomtypes ]; split it so atomtypes can go before every [ moleculetype ]
# ACPYPE atomtypes carry no atomic number; the second awk adds it (from the mass) so that
# gmx hbond (2024+) can find the ligand's N/O donors and acceptors
awk '/^\[ *atomtypes *\]/{f=1} /^\[ *moleculetype *\]/{f=0} f' LIG_GMX.itp | \
  awk '!/^ *[;[]/ && NF>=7 {m=$3; z=(m<1.5)?1:(m<13)?6:(m<15)?7:(m<17)?8:(m<20)?9:(m<31.5)?15:(m<33)?16:(m<36)?17:(m<80)?35:53; $2=$2" "z} 1' \
  > ../LIG_atomtypes.itp
awk '/^\[ *moleculetype *\]/{f=1} f' LIG_GMX.itp > ../LIG.itp
cp LIG_GMX.gro ../LIG.gro
cp posre_LIG.itp ../posre_LIG.itp      # heavy-atom restraints, 1000 kJ/mol/nm^2
cd ..

# OpenFF Sage ligand -> GROMACS files (Interchange)
from openff.toolkit import ForceField, Molecule
from openff.interchange import Interchange

lig = Molecule.from_file("lig_docked.sdf")      # coordinates = docked pose, H included
lig.name = "LIG"
sage = ForceField("openff-2.3.0.offxml")
inter = Interchange.from_smirnoff(force_field=sage, topology=[lig], box=[5, 5, 5])
inter.to_top("LIG_openff.top")
inter.to_gro("LIG_openff.gro")
# Then split it in bash; [ defaults ] is dropped, and the moleculetype part stops at [ system ]:
#   awk '/^\[ *atomtypes *\]/{f=1} /^\[ *moleculetype *\]/{f=0} f' LIG_openff.top > LIG_atomtypes.itp
#   awk '/^\[ *moleculetype *\]/{f=1} /^\[ *system *\]/{f=0} f' LIG_openff.top > LIG.itp
#   cp LIG_openff.gro LIG.gro     # atom order = lig_docked.sdf; make posre_LIG.itp with genrestr as in the CGenFF route

# CGenFF route: generate ligand heavy-atom restraints on the ligand-only coordinates
gmx editconf -f lig_gmx.pdb -o LIG.gro
printf '0 & ! a H*\nq\n' | gmx make_ndx -f LIG.gro -o index_LIG.ndx
echo 3 | gmx genrestr -f LIG.gro -n index_LIG.ndx -o posre_LIG.itp -fc 1000 1000 1000

# Append ligand atoms to the protein coordinates (same frame, nm)
python3 - <<'EOF'
p = open("protein.gro").read().splitlines()
l = open("LIG.gro").read().splitlines()
atoms = p[2:-1] + l[2:-1]
open("complex.gro", "w").write("\n".join([p[0], f"{len(atoms):5d}", *atoms, p[-1]]) + "\n")
EOF

; topol.top (AMBER + GAFF2 example; CHARMM route is identical in layout)
#include "amber14sb.ff/forcefield.itp"
#include "LIG_atomtypes.itp"          ; ligand atomtypes: after forcefield.itp, before any [ moleculetype ]

[ moleculetype ]
; Name            nrexcl
Protein_chain_A     3
...                                   ; written by pdb2gmx

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Ligand topology and its restraints (must follow the ligand [ moleculetype ])
#include "LIG.itp"
#ifdef POSRES_LIG
#include "posre_LIG.itp"
#endif

; Include water topology
#include "amber14sb.ff/tip3p.itp"
...

[ molecules ]
; Compound        #mols   (same order as complex.gro)
Protein_chain_A     1
LIG                 1

gmx editconf -f complex.gro -o box.gro -bt dodecahedron -d 1.0 -c
gmx solvate -cp box.gro -cs spc216.gro -p topol.top -o solv.gro     # OPC/TIP4P-type water: -cs tip4p.gro

gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
echo SOL | gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15

; ions.mdp  (only used to build ions.tpr for genion)
integrator      = steep
emtol           = 1000.0
emstep          = 0.01
nsteps          = 50000
nstlist         = 1
cutoff-scheme   = Verlet
coulombtype     = cutoff      ; plain cutoff: no "Ewald with net charge" warning before neutralization
rcoulomb        = 1.0
rvdw            = 1.0
pbc             = xyz

; em.mdp  (CHARMM36m settings; see the AMBER note below)
integrator      = steep
emtol           = 1000.0        ; kJ/mol/nm
emstep          = 0.01
nsteps          = 50000
nstlist         = 1
cutoff-scheme   = Verlet
coulombtype     = PME
rcoulomb        = 1.2
vdwtype         = Cut-off
vdw-modifier    = Force-switch
rvdw-switch     = 1.0
rvdw            = 1.2
DispCorr        = no
pbc             = xyz

; nvt.mdp  100 ps, protein and ligand heavy atoms restrained
define                  = -DPOSRES -DPOSRES_LIG
integrator              = md
nsteps                  = 50000         ; 50000 * 2 fs = 100 ps
dt                      = 0.002
nstxout-compressed      = 5000
nstenergy               = 500
nstlog                  = 500
continuation            = no
constraint-algorithm    = lincs
constraints             = h-bonds
lincs-iter              = 1
lincs-order             = 4
cutoff-scheme           = Verlet
nstlist                 = 20
coulombtype             = PME
rcoulomb                = 1.2
fourierspacing          = 0.16
pme-order               = 4
vdwtype                 = Cut-off
vdw-modifier            = Force-switch
rvdw-switch             = 1.0
rvdw                    = 1.2
DispCorr                = no
tcoupl                  = V-rescale
tc-grps                 = System
tau-t                   = 1.0
ref-t                   = 300
pcoupl                  = no
pbc                     = xyz
gen-vel                 = yes
gen-temp                = 300
gen-seed                = -1            ; new seed per replica

; npt.mdp  100 ps, restraints kept, C-rescale barostat
define                  = -DPOSRES -DPOSRES_LIG
integrator              = md
nsteps                  = 50000
dt                      = 0.002
nstxout-compressed      = 5000
nstenergy               = 500
nstlog                  = 500
continuation            = yes
constraint-algorithm    = lincs
constraints             = h-bonds
lincs-iter              = 1
lincs-order             = 4
cutoff-scheme           = Verlet
nstlist                 = 20
coulombtype             = PME
rcoulomb                = 1.2
fourierspacing          = 0.16
pme-order               = 4
vdwtype                 = Cut-off
vdw-modifier            = Force-switch
rvdw-switch             = 1.0
rvdw                    = 1.2
DispCorr                = no
tcoupl                  = V-rescale
tc-grps                 = System
tau-t                   = 1.0
ref-t                   = 300
pcoupl                  = C-rescale
pcoupltype              = isotropic
tau-p                   = 5.0
ref-p                   = 1.0
compressibility         = 4.5e-5
refcoord-scaling        = com           ; required with position restraints + pressure coupling
pbc                     = xyz
gen-vel                 = no

; md.mdp  100 ns production, no restraints
integrator              = md
nsteps                  = 50000000      ; 50,000,000 * 2 fs = 100 ns
dt                      = 0.002
nstxout                 = 0
nstvout                 = 0
nstfout                 = 0
nstxout-compressed      = 5000          ; one frame every 10 ps
compressed-x-grps       = System
nstenergy               = 5000
nstlog                  = 5000
continuation            = yes
constraint-algorithm    = lincs
constraints             = h-bonds
lincs-iter              = 1
lincs-order             = 4
cutoff-scheme           = Verlet
nstlist                 = 20
coulombtype             = PME
rcoulomb                = 1.2
fourierspacing          = 0.16
pme-order               = 4
vdwtype                 = Cut-off
vdw-modifier            = Force-switch
rvdw-switch             = 1.0
rvdw                    = 1.2
DispCorr                = no
tcoupl                  = V-rescale
tc-grps                 = System
tau-t                   = 1.0
ref-t                   = 300
pcoupl                  = C-rescale
pcoupltype              = isotropic
tau-p                   = 5.0
ref-p                   = 1.0
compressibility         = 4.5e-5
pbc                     = xyz
gen-vel                 = no

; AMBER14SB / AMBER19SB + GAFF2 or Sage: replace the nonbonded block with
rcoulomb                = 1.0
vdwtype                 = Cut-off
vdw-modifier            = Potential-shift
rvdw                    = 1.0
DispCorr                = EnerPres
; and delete rvdw-switch

# Optional: two coupling groups (the 2018 tutorial layout)
printf '"Protein" | "LIG"\nq\n' | gmx make_ndx -f em.gro -o index.ndx
# in nvt/npt/md.mdp:  tc-grps = Protein_LIG Water_and_ions   tau-t = 1.0 1.0   ref-t = 300 300
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr

gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt -nb gpu -pme gpu

gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
gmx mdrun -deffnm npt -nb gpu -pme gpu

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr
gmx mdrun -deffnm md -ntmpi 1 -ntomp 8 -nb gpu -pme gpu -bonded gpu -update gpu -pin on

# Three independent replicas: same minimized structure, new velocities in NVT
for r in 1 2 3; do
  mkdir -p rep$r && cd rep$r
  gmx grompp -f ../nvt.mdp -c ../em.gro -r ../em.gro -p ../topol.top -o nvt.tpr
  gmx mdrun -deffnm nvt -nb gpu -pme gpu
  gmx grompp -f ../npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p ../topol.top -o npt.tpr
  gmx mdrun -deffnm npt -nb gpu -pme gpu
  gmx grompp -f ../md.mdp -c npt.gro -t npt.cpt -p ../topol.top -o md.tpr
  gmx mdrun -deffnm md -ntmpi 1 -nb gpu -pme gpu -bonded gpu -update gpu
  cd ..
done

# 1) Interrupted run: continue from the last checkpoint (written every 15 min by default)
gmx mdrun -deffnm md -cpi md.cpt -nb gpu -pme gpu -bonded gpu -update gpu

# 2) Finished 100 ns, extend by another 100 ns (-extend is in ps)
gmx convert-tpr -s md.tpr -extend 100000 -o md_200ns.tpr
gmx mdrun -deffnm md -s md_200ns.tpr -cpi md.cpt -nb gpu -pme gpu -bonded gpu -update gpu

# Or set the end time directly: run until 200 ns
gmx convert-tpr -s md.tpr -until 200000 -o md_200ns.tpr

# 3) Queue system with a wall-time limit: stop cleanly before 23.5 h and write a checkpoint
gmx mdrun -deffnm md -cpi md.cpt -maxh 23.5

# 4) Changing mdp options while continuing (e.g. output frequency)
gmx grompp -f md_new.mdp -c npt.gro -t md.cpt -p topol.top -o md_new.tpr
gmx mdrun -s md_new.tpr -cpi md.cpt -noappend

# Make molecules whole and keep the complex centered
printf '"Protein" | "LIG"\nq\n' | gmx make_ndx -f em.gro -o index.ndx
printf 'Protein_LIG\nSystem\n' | gmx trjconv -s md.tpr -f md.xtc -n index.ndx \
    -o md_center.xtc -center -pbc mol -ur compact

# gmx_MMPBSA 1.7.x, GB with the 1.7 defaults written out explicitly
cat > mmpbsa.in <<'EOF'
&general
  sys_name="complex", startframe=5001, endframe=10000, interval=10,
  PBRadii=4,
/
&gb
  igb=8, saltcon=0.150,
/
EOF
gmx_MMPBSA -O -i mmpbsa.in -cs md.tpr -ct md_center.xtc -ci index.ndx \
    -cg Protein LIG -cp topol.top -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv

# 1) 电荷：Multiwfn 自带的脚本（Multiwfn 目录下 examples/RESP/），在 Linux 下运行
#    运行前在脚本中改好 ORCA= 与 orca_2mkl= 的路径、nprocs 和 maxcore
#    参数：结构文件 净电荷 自旋多重度；输入对接得到的带氢配体，脚本先做 B97-3c 优化再算电荷
./RESP2_ORCA.sh lig.mol2 0 1
#    结构已经优化过时，用不做优化的版本
./RESP2_ORCA_noopt.sh lig.mol2 0 1
#    有 Gaussian 时用 RESP2.sh 或 RESP.sh，用法相同
#    输出 lig.chg：前四列是元素和坐标，最后一列是电荷，原子顺序与输入文件一致

# 2) 拓扑：启动 ./sobtop 后依次输入（Sobtop 主页例 2 的 GAFF 路线）
#    lig.mol2        对接构象、带氢、键级正确的 mol2
#    lig.chg         在主菜单直接输入 chg 文件路径，载入 RESP2 电荷
#    2, [回车]       写出 gro，坐标取自 mol2
#    1, 2, 4         写 GROMACS 拓扑；指认 GAFF 原子类型；从参数库取参数，缺失项由程序猜
#    [回车], [回车]  使用默认的 top 和 itp 路径