摘要本文系统梳理 Amber 分子动力学模拟蛋白-核酸复合物蛋白质-DNA / 蛋白质-RNA的完整流程从 PDB 结构与 pdb4amber 预处理入手重点说明 ff19SB OL15/OL21/OL3 OPC JC 离子的力场组合随后演示 tleap 两段式 combine 构建体系的完整命令接着给出强约束最小化、分段升温、弱约束平衡与生产模拟的可复用输入文件并总结 pmemd 运行时的常见坑最后以 cpptraj 实战为验收核心介绍 nastruct 螺旋参数、糖环 pucker 等核酸特征分析方法补充离子浓度、Mg²⁺、周期性边界等特殊注意事项及常见错误排查。全文以 Drew-Dickerson 十二聚体为官方教程基线命令参数可直接照抄复用。全文约 10 千字 预计阅读 25 分钟本文介绍Amber分子动力学模拟结构准备蛋白-核酸复合物蛋白质–DNA / 蛋白质–RNA的分子动力学模拟操作及注意事项。文中命令与输入文件参数取自 Amber 官方核酸模拟教程AMBER HubDrew-Dickerson 十二聚体案例与近年力场评测文献可直接照抄复用。做转录因子–DNA、Cas9–sgRNA–DNA、核酸酶这类体系的模拟时很多人直接套用纯蛋白的流程结果跑出来 DNA 解旋、能量爆炸或者 RMSD 永不收敛。蛋白–核酸复合物与纯蛋白体系有三个本质差异每一步操作都要围绕它们展开强负电体系每两个核苷酸就有一个磷酸基团一条 12 bp 双链 DNA 净电荷 -22离子环境直接决定结构稳定性力场组合蛋白力场 核酸力场必须同时加载且核酸力场的版本选择比蛋白更讲究选错会系统性破坏螺旋几何构象验证DNA 有 A/B/Z 三种构象RNA 以 A-form 为主糖环 puckerC2′-endo vs C3′-endo是否正确是判断模拟可信度的硬指标纯蛋白流程里没有这一步。相关教程与核心文献官方教程教程内容与本文关系AMBER 官方教程总目录全部基础/高级教程索引按体系类型检索后续进阶教程的入口AMBER Hub 核酸模拟教程Drew-Dickerson 十二聚体全流程本文 tleap/mdin/nastruct 命令的直接出处核心文献文献为什么值得先读Tian et al., ff19SB, JCTC 2019. DOI 10.1021/acs.jctc.9b00591本文蛋白力场选择的原始依据——理解 CMAP 项为何必须配 OPC 水Izadi Onufriev, OPC 水模型, JPCL 2014. DOI 10.1021/jz501780a核酸静电描述更准的水模型与 JC 离子参数配套加载的原理OL15 力场系列综述, JCTC 2016. DOI 10.1021/acs.jctc.6b00186弄清 ff99 → bsc0 → OL15 的骨架二面角修正脉络避免误用旧参数Love et al., JCTC 2023. DOI 10.1021/acs.jctc.3c01164提出 dsDNA 用 OL21 OPC 的新推荐决定你该加载哪个 leaprc一、准备工作1. 获取初始结构来源PDBX 射线或 Cryo-EM 结构例如1ABC.pdb。检查清单是否包含完整蛋白和核酸链是否有缺失残基/原子可用Modeller、ChimeraX补全非标准残基修饰碱基、甲基化碱基、甲基化氨基酸需单独参数化参见非标准残基处理指南删除无关水分子、配体除非必要。✅ 建议保留结合位点附近的结晶水介导蛋白–核酸识别的水桥但后续统一用显式水模型重建溶剂环境。2. pdb4amber 预处理pdb4amber protein_dna.pdb --dry --reduce --overwrite--dry移除水可选--reduce修复 HIS 命名HID/HIE并补氢自动标准化残基命名。⚠️实测坑部分 PDB 条目里核酸链整链是HETATM记录而不是ATOM。--dry只按记录类型删水但如果后续手动按HETATM清理杂原子会把整条 DNA 一起删掉。预处理后务必grep -c DA\|DT\|DG\|DC *.pdb确认核酸链还在。二、选择合适的力场核酸模拟的地基Amber 官方推荐组合基于 Amber20/22 推荐实践与 2023 年力场评测组分力场加载命令tleap说明蛋白质ff19SBsource leaprc.protein.ff19SB配 OPC 水模型标定DNAOL15教程默认或OL212023 评测推荐source leaprc.DNA.OL15bsc1 为备选RNAOL3source leaprc.RNA.OL3RNA 标准选择水模型OPC推荐或 TIP3Psource leaprc.water.opcOPC 对核酸静电描述更准离子JC ions匹配所用水模型自动加载随 leaprc.water.* 加载不要跨水模型混用离子参数三个要点OL15 是什么。加载leaprc.DNA.OL15时终端会打印一行说明OL15 force field for DNA (99bsc0-betaOL1-eps-zetaOL1-chiOL4)——即在 ff99 基础上叠加 bsc0 的 α/γ 修正、βOL1、ε/ζOL1、χOL4 四组骨架二面角修正。这些修正专门解决早期力场把 B-DNA 慢慢扭曲成非物理构象的问题。2023 年之后的最佳实践。Love 等J. Chem. Theory Comput. 2023对 Amber 核酸力场做系统 vdW 参数扫描后推荐 dsDNA 用 OL21 OPC 水模型。AmberTools 21 起提供leaprc.DNA.OL21用 AmberTools 20 或更早版本则用 OL15。蛋白与水模型要配套。ff19SB 的 CMAP 项是配合 OPC 水优化标定的官方教程包括核酸教程现在统一source leaprc.water.opcOPC 与 JC 离子参数随之自动加载。TIP3P 仍可用传统组合 ff14SBTIP3PJC但不要让蛋白用 ff19SBOPC、离子却手动塞 TIP3P 的参数。⚠️ 不要混用旧力场如 ff99 parmbsc0 的原始参数可能导致能量不一致与螺旋畸变。三、使用 tleap 构建体系3.1 纯核酸热身Drew-Dickerson 十二聚体官方教程实录Amber 官方核酸教程用 d(CGCGAATTCGCG)Drew-Dickerson dodecamerDDD做示范——这是核酸模拟领域被研究最多的短双链从 A-DNA 构象出发观察模拟中自发转变为 B-DNA。完整 tleap 会话$ tleap source leaprc.DNA.OL15 source leaprc.water.opc dna loadpdb A-dna.pdb Loading PDB file: ./A-dna.pdb total atoms in file: 758 check dna Checking dna.... WARNING: The unperturbed charge of the unit: -22.000000 is not zero. check: Warnings: 1 Unit is OK. addions dna Na 0 solvateoct dna OPCBOX 9.0 saveamberparm dna dna.prmtop dna.rst7 quit逐行解读这些数字就是你的对照标准758 个原子12 bp 双链 DNA 的重原子 末端残基标准处理后的规模净电荷 -2222 个磷酸基团check报这个 WARNING 是正常预期不是错误——看到它说明磷酸骨架完整如果电荷不是 -2×(链长-1) 的量级说明缺了核苷酸addions dna Na 0加 Na⁺ 直到净电荷为 0中和离子优先放在负电位点附近solvateoct dna OPCBOX 9.0 truncated octahedron 八面体盒子溶质到盒边至少 9 Å教程体系约加 5336 个水分子。八面体比立方盒省相当数量的水省计算量长核酸体系推荐saveamberparm输出拓扑dna.prmtop 坐标dna.rst7。3.2 蛋白 B-DNA 复合物构建两段式 combine蛋白–核酸复合物推荐把蛋白链与核酸链拆成独立 PDB、分别加载再 combineGENESIS 教程用 PDB 3LNQ 演示的同一套思路Amber 下通用# 加载力场蛋白 DNA 水三行都要 source leaprc.protein.ff19SB source leaprc.DNA.OL15 source leaprc.water.opc # 若含 RNAsource leaprc.RNA.OL3 # 蛋白与核酸分别加载残基命名须符合 Amber 规范 # 标准 DNA 残基名DA5/DA/DA3, DG, DC, DT5/3 端带编号 # 标准 RNA 残基名A5/A/A3, G, C, U protein loadpdb protein_clean.pdb dna loadpdb dna_clean.pdb complex combine { protein dna } # 检查警告重点关注 unknown residue 或 missing atoms check complex desc complex # 中和净电荷核酸体系通常加大量 Na addIons complex Na 0 # 追加生理盐浓度按盒子体积估算离子对数如 150 mM NaCl addIons complex Na 12 addIons complex Cl- 12 # 溶剂化最小水层 10–12 Å solvateoct complex OPCBOX 10.0 # 保存拓扑和坐标 saveamberparm complex system.prmtop system.inpcrd savepdb complex system_solvated.pdb quit关键检查点若 PDB 中 DNA 残基名为ADE,THY等旧命名需重命名为DA,DT否则 tleap 报 unknown residuecombine之后 tleap 会把两条链重排成连续编号蛋白 1..N核酸 N1..M后面写 restraintmask 时以重排后的编号为准用desc complex核对5′/3′ 末端残基类型DA5/DA3没给对时 tleap 会在末端缺磷酸或多余氧check阶段就能暴露。四、能量最小化与逐步平衡官方教程对核酸体系用的是强约束最小化 → 约束下缓慢升温 → 弱约束 NPT 平衡三步比纯蛋白流程更保守核酸初始结构常有实验几何与力场不匹配的张力。三个输入文件可直接复用mini.in —— 25 kcal/mol 强约束下最小化 DNAMinimzing the system with 25 kcal/mol restraints on DNA, 500 steps of steepest descent and 500 of conjugated gradient cntrl imin1, ntx1, irest0, ntpr50, ntf1, ntb1, cut9.0, ntr1, maxcyc1000, ncyc500, ntmin1, restraintmask:1-24, restraint_wt25.0, end ewald ew_type 0, skinnb 1.0, endheat.in —— 约束下 100 K → 300 K 分段升温NVTHeating the system with 25 kcal/mol restraints on DNA, Vconst cntrl imin0, ntx1, ntpr500, ntwr500, ntwx500, ntwe500, nscm5000, ntf2, ntc2, ntb1, ntp0, nstlim50000, t0.0, dt0.002, cut9.0, tempi100.0, ntt1, ntr1, nmropt1, restraintmask:1-24, restraint_wt25.0, end wt typeTEMP0, istep10, istep25000, value1100.0, value2300.0, end wt typeTEMP0, istep15001, istep250000, value1300.0, value2300.0, end wt typeEND, endeq.in —— 0.5 kcal/mol 弱约束下 NPT 平衡 500 psEquilibrating the system with 0.5 kcal/mol restraints on DNA, during 500 ps, at constant T300 K, P1 ATM cntrl imin0, ntx5, ntpr500, ntwr500, ntwx500, ntwe500, nscm5000, ntf2, ntc2, ntb2, ntp1, tautp0.2, taup0.2, nstlim25000, t0.0, dt0.002, cut9.0, ntt1, ntr1, irest1, restraintmask:1-24, restraint_wt0.5, end ewald ew_type 0, skinnb 1.0, end三个文件的要点约束从 25 → 0.5 kcal/mol 阶梯释放最小化和升温阶段把 DNA:1-24为教程体系 24 个核苷酸残基位置约束住只让水和离子松弛避免溶剂化冲击破坏晶体几何平衡阶段降到 0.5 让核酸开始弛豫但不上蹿下跳。生产前再把ntr0放开跑一段无约束平衡升温用 nmropt1 TEMP0 分段前 5000 步从 100 K 线性升到 300 K之后恒温——比一步到位加热温和得多ntc2/ntf2SHAKE从升温阶段就启用配合 dt0.002 ps。运行顺序注意每步的-ref都要给pmemd -O -i mini.in -o mini.out -p dna.prmtop -c dna.rst7 -r mini.rst -ref dna.rst7 pmemd -O -i heat.in -o heat.out -p dna.prmtop -c mini.rst -r heat.rst -x heat.nc -e heat.mden -ref mini.rst pmemd -O -i eq.in -o eq.out -p dna.prmtop -c heat.rst -r eq.rst -x eq.nc -e eq.mden -ref heat.rst⚠️实测坑pmemd 16 运行时一律加-Ooverwrite否则输出文件已存在时直接报Unit NN Error on OPEN终止设置了ntwe的输入文件必须配-e mden文件名否则同样报 OPEN 错误。教程网页上的命令行没写这两点照抄会翻车。蛋白–核酸复合物在 eq 之后建议再加一步完全无约束的短平衡ntr050–100 ps观察蛋白–核酸界面接触是否维持再进生产。五、生产模拟官方教程的 md.inNPT200 ps 量级示例production cntrl imin0, ntx5, ntpr1000, ntwr1000, ntwx1000, ntwe1000, nscm1000, ntf2, ntc2, ntb2, ntp1, tautp5.0, taup5.0, nstlim100000, t0.0, dt0.002, cut9.0, ntt1, irest1, iwrap1, ioutfm1, end ewald ew_type 0, skinnb 1.0, endpmemd -i md.in -o md.out -p dna.prmtop -c eq.rst -r md.rst7 -x md.nc -inf mdinfo生产实践建议教程的nstlim100000 × dt0.002 200 ps只是演示量级蛋白–核酸复合物的科研模拟建议 100 ns 起步蛋白–核酸界面水合与离子再分布的弛豫时间比蛋白–蛋白界面长大体系上 GPUpmemd.cuda直接替换pmemd输入文件不用改tautp/taup5.0比平衡阶段0.2放松避免过度耦合扭曲动力学ioutfm1输出 NetCDF 轨迹长模拟必须ASCII mdcrd 会大到没法处理。六、核酸特征分析cpptraj 实战这是蛋白–核酸流程区别于纯蛋白流程的核心验收环节。官方教程用 cpptraj 完成三类分析6.1 轨迹预处理parm dna.prmtop trajin md.nc autoimage strip :WAT trajout md-nowater.nc netcdfautoimage处理周期性边界strip :WAT去水减小轨迹体积。去水后要用 parmed 生成对应的新拓扑strip :WATparmout dna-nowater.prmtop才能在 VMD 里看。6.2 螺旋参数nastructcpptraj 的nastruct命令一次输出核酸构象的三类参数文件输出文件内容列BP.data碱基对参数Shear / Stretch / Stagger / Buckle / Propeller / Opening / HB / Major / MinorBPstep.data碱基步参数Shift / Slide / Rise / Tilt / Roll /TwistHelix.data螺旋轴参数X-disp / Y-disp / Rise / Incl. / Tip / Twist判断标准B-DNA 常见值域Rise ≈ 3.3–3.4 ÅTwist ≈ 34–36°Slide 为小幅负值。若 Twist 漂到 30° 以下或 Rise 超过 3.6 Å通常意味着力场/离子环境出了问题。教程里用一段 awk 提取中心碱基步的 Twist 随时间演化sed -e /^$/d BPstep.data_6-7 BPstep.data_6-7.new awk {print $1 \t $9} BPstep.data_6-7.new twist_6-7.datDDD 案例的示范结果A-DNA 起点在模拟中自发向 B-DNA 转变Twist 在约第 180 帧附近抬升与 RMSD 收敛时间点对应——构象转变事件要同时用 RMSD 螺旋参数交叉确认只看 RMSD 会漏判。6.3 糖环 puckercpptraj 的pucker命令计算五元糖环的伪旋转相位角Altona Sundaralingam 方法。注意 mask 写法一条:1-24C1 ...的范围 mask 会被 cpptraj 静默截成只算第一个残基每列 mask 只取首个匹配原子必须逐残基写pucker p1 :1C1 :1C2 :1C3 :1C4 :1O4 out pucker.dat pucker p2 :2C1 :2C2 :2C3 :2C4 :2O4 out pucker.dat ...逐残基写全C2′-endo 对应伪旋转相位角约 110–180°B-DNA 特征C3′-endo 约 0–40°A-DNA / RNA 特征。实测示例Drew-Dickerson 十二聚体A-DNA 起点OL15OPC按本节流程跑 200 ps 生产全链 24 个残基 × 100 帧中仅 31% 采样落在 C2′-endo 区残基 2–6 和 15–19两条链的 5′ 端半段超过 80% 的帧仍停留在 C3′-endo 区——A→B 转变在 200 ps 内只完成了一半。这与官方教程~800 帧太少、延长模拟分布才平滑的提示一致糖环 pucker 的构象转变是慢变量判断 B-form 是否达成需要 ns 级采样200 ps 只够看到转变开始。RNA 体系判据反过来主峰应落在 C3′-endo。6.4 其他常规分析任务工具轨迹分析cpptraj,pytraj氢键/接触cpptraj nativecontacts / hbond蛋白–DNA 界面 RMSDcpptraj rms分链分别算别用全体系一把算结合自由能MMPBSA.py, TI/FEP精细螺旋参数Curves,3DNA外部专用工具cpptraj 交叉验证用七、特殊注意事项1. 离子浓度设置蛋白–核酸复合物带强负电磷酸骨架需足够阳离子屏蔽。先addIons ... 0中和再按盒子体积追加 100–150 mM NaCl高盐体系的离子数要自己算n ≈ 0.15 × 0.6022 × V(ų) 个离子对别拍脑袋可用parmed查看残基电荷核对净电荷parmed system.prmtop parmed printDetails :1-1002. Mg²⁺ 处理若结构中有 Mg²⁺常见于核酸酶、核酶、聚合酶Amber 的 JC 离子模型支持 Mg²⁺但对配位键的描述精度有限——12-6 LJ 球形模型无法重现八面体配位几何对关键结构/催化 Mg²⁺建议用Li/Merz 12-6-4 非键模型在 12-6 基础上加诱导偶极项配位几何明显改善或MCPB.py键合模型成键参数化最准但最贵弱结合的游离 Mg²⁺只是电荷屏蔽用用 JC 模型即可。3. 周期性边界与盒子大小核酸较长且各向异性长条形旋转扩散快需确保任意方向水层 ≥10 Å避免双链首尾穿过周期边界与自身相互作用solvateOct八面体盒子对棒状溶质特别省水——教程 DDD 体系 9 Å 水层只用约 5336 个水立方盒要多出一截。4. 验证核酸构象验收清单平衡与生产后逐项检查 DNA/RNA 是否保持目标构象B/A/Z 或 A-formnastruct螺旋参数落在目标构象值域见 6.2糖环 pucker 主峰位置正确见 6.3RMSD 收敛后检查是否发生 A↔B 构象漂移DDD 是有意让它转你的体系如果不是研究对象就该警惕。八、常见错误排查问题解决方案tleap 报 Unknown residue XXX检查残基命名ADE→DA 等用pdb4amber标准化5′/3′ 末端残基名DA5/DA3是否给对check报净电荷异常教程 DDD 报 -22 是正常22 个磷酸基团电荷偏离 -2×(残基数-1) 才说明缺链/缺残基能量爆炸minimization 失败检查原子重叠提高第一阶段约束强度25 kcal/mol 起步确认没有双占据的晶体构象副本DNA 解旋/断裂/末端 fray力场是否为 OL15/OL21/OL3是否缺中和离子水层是否够厚导致自相互作用Twist/Rise 漂出 B-DNA 值域检查是否误用 ff99 旧参数换 OPC 水模型核对离子参数与水模型匹配RMSD 永不收敛先看核酸部分单独的 RMSD——若只有核酸漂多半是构象转变nastruct 确认若整体漂检查盒子/温度耦合pdb4amber 后核酸链消失核酸链是 HETATM 记录被当杂原子清理重新处理按链名保留九、结语蛋白–核酸复合物模拟的完整心法一句话力场组合蛋白 ff19SB 核酸 OL15/OL21/OL3 OPC 水 匹配离子→ tleap 两段式 combine → 强约束起步的阶梯式平衡 → 用 nastruct 和 pucker 验收构象。前三点决定模拟能不能稳定跑下去最后一点决定结果能不能信。如提供具体体系如 p53–DNA 复合物、Cas9–sgRNA–DNA 三元复合物可进一步定制流程如处理 sgRNA 修饰、多链组装、Mg²⁺ 配位等。关键字蛋白-核酸复合物、分子动力学模拟、AMBER、tleap、OL15力场、OPC水模型、nastruct、糖环pucker参考文献与延伸阅读AMBER Hub 官方核酸模拟教程Drew-Dickerson dodecamer 全流程本文 tleap/mdin/nastruct 命令来源: https://amberhub.chpc.utah.edu/analisis-of-nucleic-acid-simulation/OL15 力场系列综述之一J. Chem. Theory Comput., 2016: DOI 10.1021/acs.jctc.6b00186核酸力场变体系统评估Galindo-Murillo et al., Nucleic Acids Res. 2017, 45(7): 4217: DOI 10.1093/nar/gkx136Love et al. dsDNA 力场 vdW 参数扫描提出 OL21OPC 推荐J. Chem. Theory Comput., 2023: DOI 10.1021/acs.jctc.3c01164OPC 水模型Izadi, Anandakrishnan, Onufriev, J. Phys. Chem. Lett. 2014, 5, 3863: DOI 10.1021/jz501780aOL15/OL21 官方主页Olomouc 力场组: http://ffol.upol.cz