Amber99SB-ILDN力场MD模拟mdp文件及数据处理脚本分享
在我的文章《在云服务器AutoDL实现分子动力学模拟全流程》中我分享了MD的步骤和相关的命令行，而本文中我将分享其中提到的mdp文件和python绘图脚本。这一部分涉及非常多可选参数，我将进行注释。主要由AI生成，这一部分涉及的知识太多，先应用为主，只在一些适应性修改上做讨论
ions.mdp
在VS code中新建文档，复制进去保存.mdp格式即可，后面同理
integrator = steep nsteps = 0 cutoff-scheme = Verlet constraint_algorithm = LINCSminim.mdp
; 能量最小化参数 —— 适用于蛋白质-蛋白质复合体 + 离子体系 integrator = steep ; 算法：最陡下降法（初始结构松弛最稳定） emtol = 1000.0 ; 收敛条件：最大作用力低于 1000 kJ/mol/nm 就停止 emstep = 0.01 ; 能量最小化步长（默认值，安全可靠） nsteps = 50000 ; 最大步数，防止不收敛时无限运行 ; 近邻搜索与边界条件 cutoff-scheme = Verlet ; 近邻搜索方案 nstlist = 10 ; 每10步更新一次近邻列表 ns_type = grid ; 用网格方式搜索近邻，速度快、稳定 pbc = xyz ; 开启三维周期性边界条件（溶液模拟必须开） ; 静电相互作用与范德华作用 coulombtype = PME ; 粒子网格Ewald，处理离子、盐溶液必需，更准确 rcoulomb = 1.2 ; 静电作用截断距离（1.2 适合 AMBER 力场） rvdw = 1.2 ; 范德华作用截断距离（与静电保持一致） DispCorr = EnerPres ; 长程色散校正，对蛋白-蛋白结合、疏水作用更准确 ; 约束设置 constraints = none ; 初始能量最小化不加任何约束，让结构完全松弛nvt.mdp
; nvt.mdp —— 作为grompp的输入文件，生成动力学运行文件nvt.tpr ; 参数适配卵母细胞模拟环境下的蛋白-蛋白复合物体系 ; ====================== 预处理参数 ====================== define = -DPOSRES ; 开启蛋白位置约束（固定蛋白骨架） ; ====================== 模拟运行控制 ====================== integrator = md ; 积分算法：蛙跳算法Leap-frog分子动力学 nsteps = 50000 ; 总模拟步数50000步；2 fs/步，总时长100000 fs = 100 ps dt = 0.002 ; 积分时间步长，单位纳秒，即2飞秒(2 fs) ; ====================== 输出文件打印控制 ====================== nstxout = 500 ; 每500步输出坐标轨迹，间隔1.0 ps nstvout = 500 ; 每500步输出原子速度，间隔1.0 ps nstenergy = 500 ; 每500步输出能量数据，间隔1.0 ps nstlog = 500 ; 每500步更新日志log文件，间隔1.0 ps ; ====================== 邻域搜索 & 短程相互作用 ====================== cutoff-scheme = Verlet ; 邻域截断方案：Verlet缓冲邻域列表（GPU高效） ns_type = grid ; 邻域搜索方式：网格划分搜索 nstlist = 10 ; 每10步更新一次邻域列表（20 fs更新一次） rcoulomb = 1.2 ; 静电短程相互作用截断半径，单位nm rvdw = 1.2 ; 范德华短程相互作用截断半径，单位nm ; ====================== 长程静电计算 ====================== coulombtype = PME ; 长程静电算法：粒子网格埃瓦尔德PME pme_order = 4 ; PME插值阶数，4阶立方插值 fourierspacing = 0.16 ; FFT傅里叶网格间距，单位nm ; ====================== 温度耦合（控温） ====================== tcoupl = V-rescale ; 控温算法：改良Berendsen速度标度恒温器 tc-grps = Protein Non-Protein ; 分两组控温：蛋白、非蛋白（溶剂/配体等） tau_t = 0.1 0.1 ; 两组控温弛豫时间常数，单位ps ref_t = 310 310 ; 两组目标参考温度，310 K（人体生理温度） ; ====================== 压力耦合（控压） ====================== pcoupl = no ; NVT系综不开启压力耦合，体积固定 ; ====================== 周期性边界条件 ====================== pbc = xyz ; X/Y/Z三维均开启周期性边界 ; ====================== 范德华截断色散校正 ====================== DispCorr = EnerPres ; 对截断范德华作用同时校正能量与压力 ; ====================== 初始速度生成 ====================== gen_vel = yes ; 初始化原子速度，按麦克斯韦速率分布分配 gen_temp = 310 ; 生成初速度对应的参考温度310 K gen_seed = -1 ; 随机种子=-1，每次运行自动生成随机种子npt.mdp
npt是恒压步，顾名思义容易出现体系崩坏，我注释了三处可修改选项，如果还是不行，可以从pdb开始把两个模型拉开几埃，重新把各步骤做一遍。
; npt.mdp —— 微管蛋白复合物压力平衡模拟输入文件 ; ====================== 预处理参数 ====================== ; 如果出现了原子间作用过强的报错，请用分号注释掉define这行，并启用下面额外注释的修改 define = -DPOSRES ; 对蛋白施加位置约束，限制蛋白整体大幅位移 ; ====================== 模拟运行控制 ====================== integrator = md ; 积分器：蛙跳Leap-frog分子动力学算法 nsteps = 50000 ; 总模拟步数50000步，总时长100 ps dt = 0.002 ; 积分步长0.002 ns，即2 fs ; ====================== 输出打印控制 ====================== nstxout = 500 ; 每500步输出坐标轨迹，间隔1.0 ps nstvout = 500 ; 每500步输出原子速度，间隔1.0 ps nstenergy = 500 ; 每500步输出能量文件，间隔1.0 ps nstlog = 500 ; 每500步更新运行日志，间隔1.0 ps ; ====================== 邻域搜索与截断半径（CHARMM36m力场标准参数） ====================== cutoff-scheme = Verlet ; Verlet缓冲邻域列表方案，GPU计算效率更高 ns_type = grid ; 网格法邻域原子搜索 nstlist = 10 ; 每10步更新一次邻域列表 rcoulomb = 1.2 ; 静电短程相互作用截断半径 1.2 nm rvdw = 1.2 ; 范德华短程相互作用截断半径 1.2 nm ; ====================== 长程静电计算 ====================== coulombtype = PME ; 粒子网格埃瓦尔德算法，高适配高氯化钾离子环境 pme_order = 4 ; PME 4阶立方插值 fourierspacing = 0.16 ; FFT傅里叶网格间距 0.16 nm ; ====================== 温度耦合（控温） ====================== tcoupl = V-rescale ; V-rescale速度标度恒温器 tc-grps = Protein Non-Protein ; 分两组控温：蛋白、非蛋白（溶剂/离子等） tau_t = 0.1 0.1 ; 两组温度弛豫时间常数，单位ps ref_t = 310 310 ; 目标参考温度310 K ; ====================== 压力耦合（NPT新增核心参数） ====================== ; pcoupl可换为berdensen，更温和 pcoupl = C-rescale ; 现代稳定型恒压算法，C-rescale压力标度 pcoupltype = isotropic ; 各向同性控压：三维盒子同步缩放 ; tau_p可减为2.0减缓升压速度，防止overlap的原子弹开 tau_p = 5.0 ; 压力弛豫时间常数，单位ps ref_p = 1.0 ; 目标参考压力 1 bar（标准大气压） compressibility = 4.5e-5 ; 水的等温压缩系数 ; ====================== 周期性边界条件 ====================== pbc = xyz ; X/Y/Z三维均开启周期性边界 ; ====================== 范德华色散校正 ====================== DispCorr = EnerPres ; 对截断范德华作用同时校正能量与压力 ; ====================== 初始速度与续跑设置 ====================== gen_vel = no ; 不重新生成初始速度，读取NVT平衡轨迹中的速度 continuation = yes ; 续跑模式，接续上一步NVT模拟结果运行 ; ====================== 键约束（NPT新增参数） ====================== constraints = h-bonds ; 约束所有含氢原子的化学键长度，允许更大步长 constraint_algorithm = lincs ; LINCS约束算法（GROMACS标准蛋白约束算法） lincs_iter = 1 ; LINCS迭代修正次数 lincs_order = 4 ; LINCS高阶展开阶数md.mdp
; md.mdp —— 用于微管蛋白二聚体结合相互作用分析的生产分子动力学模拟 ; ====================== 预处理参数 ====================== ; NO define = -DPOSRES. ; ====================== 模拟运行控制 ====================== integrator = md ; 积分算法：蛙跳Leap-frog分子动力学 nsteps = 50000000 ; 总步数50000000步，总模拟时长100 ns dt = 0.002 ; 积分时间步长0.002 ns = 2 fs ; ====================== 输出打印控制（生产模拟精简输出，节省硬盘） ====================== nstxout = 0 ; 不输出高精度未压缩坐标文件，节省存储空间 nstvout = 0 ; 不输出原子速度文件 nstfout = 0 ; 不输出原子受力文件 nstxout-compressed = 5000 ; 每5000步输出压缩轨迹，间隔10.0 ps compressed-x-grps = System ; 压缩轨迹输出整个体系所有原子 nstenergy = 5000 ; 每5000步输出能量数据，间隔10.0 ps nstlog = 5000 ; 每5000步更新运行日志文件，间隔10.0 ps ; ====================== 邻域搜索与相互作用截断参数 ====================== cutoff-scheme = Verlet ; Verlet缓冲邻域列表，GPU高效计算 ns_type = grid ; 网格划分法搜索邻近原子 nstlist = 10 ; 每10步刷新一次邻域原子列表 rcoulomb = 1.2 ; 静电短程相互作用截断半径 1.2 nm rvdw = 1.2 ; 范德华短程相互作用截断半径 1.2 nm ; ====================== 长程静电计算 ====================== coulombtype = PME ; 粒子网格埃瓦尔德算法，适配水溶液离子体系 pme_order = 4 ; PME 4阶立方插值 fourierspacing = 0.16 ; FFT傅里叶网格间距 0.16 nm ; ====================== 温度耦合控温 ====================== tcoupl = V-rescale ; V-rescale速度标度恒温器 tc-grps = Protein Non-Protein ; 分两组控温：蛋白、非蛋白（溶剂/离子） tau_t = 0.1 0.1 ; 两组温度弛豫常数，单位ps ref_t = 310 310 ; 目标模拟温度 310 K ; ====================== 压力耦合（切换为Parrinello-Rahman） ====================== pcoupl = Parrinello-Rahman ; 适合精准热力学性质计算的恒压算法 pcoupltype = isotropic ; 各向同性控压，三维模拟盒同步缩放 tau_p = 5.0 ; 压力弛豫时间常数 5.0 ps ref_p = 1.0 ; 目标参考压力 1 bar（标准大气压） compressibility = 4.5e-5 ; 液态水的等温压缩系数 ; ====================== 周期性边界条件 ====================== pbc = xyz ; X/Y/Z三个方向全部开启周期性边界 ; ====================== 范德华截断色散校正 ====================== DispCorr = EnerPres ; 同时校正能量与压力的范德华色散校正 ; ====================== 初始速度与续跑设置 ====================== gen_vel = no ; 不重新随机生成初速度，读取NPT平衡轨迹中的速度 continuation = yes ; 续跑模式，接续NPT平衡后的体系继续模拟 ; ====================== 化学键约束 ====================== constraints = h-bonds ; 约束全部含氢化学键长度，稳定使用2 fs步长 constraint_algorithm = lincs ; LINCS约束算法（蛋白体系标准约束方案） lincs_iter = 1 ; LINCS约束修正迭代次数 lincs_order = 4 ; LINCS高阶展开阶数contacts计算
VS code存成.py文件
import numpy as np import matplotlib.pyplot as plt # 1. 加载数据 # GROMACS 的 xvg 文件包含以 # 或 @ 开头的注释行 try: data = np.loadtxt("num_contacts.xvg", comments=["#", "@"]) time = data[:, 0] / 1000 # 假设原始单位是 ps，转成 ns contacts = data[:, 1] except Exception as e: print(f"读取文件出错: {e}") exit() # 2. 绘图美化 plt.figure(figsize=(8, 5), dpi=100) plt.plot(time, contacts, color='#2E7D32', linewidth=1.5, label='Atomic Contacts') # 3. 添加平滑曲线（可选，用于观察趋势） #if len(contacts) > 50: # window = int(len(contacts) / 20) # smoothed = np.convolve(contacts, np.ones(window)/window, mode='same') # plt.plot(time, smoothed, color='#FF5722', linestyle='--', label='Trend (Moving Avg)') # 4. 图表细节设定 plt.xlabel("Time (ns)", fontsize=12) plt.ylabel("Number of Contacts", fontsize=12) plt.title("A and B Protein Interaction Dynamics", fontsize=14) plt.legend(loc='best') plt.grid(True, linestyle=':', alpha=0.6) plt.savefig("contact_ab.png", dpi=300, bbox_inches='tight') # 自动调整布局防止标签遮挡 plt.tight_layout() plt.show()Hbonds计算
import numpy as np import matplotlib.pyplot as plt # 读取氢键数据 data = np.loadtxt("hbonds_interface.xvg", comments=["@", "#"]) time_ps = data[:, 0] time_ns = time_ps / 1000 # 转成纳秒 hbond_num = data[:, 1] # 画图 plt.figure(figsize=(10, 4)) plt.plot(time_ns, hbond_num, color='darkred', linewidth=1.2) plt.xlabel("Time (ns)") plt.ylabel("Number of interchain H-bonds") plt.title("Hydrogen bonds between Chain A and Chain B") plt.grid(alpha=0.3) plt.ylim(bottom=0) plt.savefig("hbond_ab.png", dpi=300, bbox_inches='tight') plt.show() # 打印关键统计 print("平均氢键数:", np.mean(hbond_num)) print("最小氢键数:", np.min(hbond_num)) print("最大氢键数:", np.max(hbond_num))