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.mdpnpt是恒压步顾名思义容易出现体系崩坏我注释了三处可修改选项如果还是不行可以从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), dpi100) plt.plot(time, contacts, color#2E7D32, linewidth1.5, labelAtomic Contacts) # 3. 添加平滑曲线可选用于观察趋势 #if len(contacts) 50: # window int(len(contacts) / 20) # smoothed np.convolve(contacts, np.ones(window)/window, modesame) # plt.plot(time, smoothed, color#FF5722, linestyle--, labelTrend (Moving Avg)) # 4. 图表细节设定 plt.xlabel(Time (ns), fontsize12) plt.ylabel(Number of Contacts, fontsize12) plt.title(A and B Protein Interaction Dynamics, fontsize14) plt.legend(locbest) plt.grid(True, linestyle:, alpha0.6) plt.savefig(contact_ab.png, dpi300, bbox_inchestight) # 自动调整布局防止标签遮挡 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, colordarkred, linewidth1.2) plt.xlabel(Time (ns)) plt.ylabel(Number of interchain H-bonds) plt.title(Hydrogen bonds between Chain A and Chain B) plt.grid(alpha0.3) plt.ylim(bottom0) plt.savefig(hbond_ab.png, dpi300, bbox_inchestight) plt.show() # 打印关键统计 print(平均氢键数:, np.mean(hbond_num)) print(最小氢键数:, np.min(hbond_num)) print(最大氢键数:, np.max(hbond_num))