简介本资源是一套面向计算化学、材料科学领域科研人员与高年级研究生的批量分子动力学MD模拟自动化工具集旨在解决重复性MD任务设置繁琐、参数管理低效、多条件对比分析困难等实际问题。压缩包共427个文件14.62MB涵盖130个AMBER输入控制文件.in、37个分子结构文件.pdb、24组力场参数体系含.prmtop、.frcmod、.lib、.inpcrd、.mol2等以及17个Shell脚本.sh和1个核心Python脚本.py完整支撑从分子预处理如ANTECHAMBER电荷生成、系统构建、GROMACS/LAMMPS兼容参数转换到批量提交与日志解析的全流程。已有694人学习下载资源以virtual_MD为核心目录组织结构清晰、示例即用附带典型log输出与inf状态监控文件便于快速上手、调试验证及拓展定制是开展温度/压力扫描、构象采样或溶剂化效应研究的高效实践基线。1. 批量分子动力学模拟脚本和示例不是“一键跑完”而是把 200 个体系从建模→提交→解析结果全链路串起来的工程化工具包你有没有试过手头有 37 个配体-蛋白复合物每个要扫 5 个温度点、3 种离子浓度光是写 LAMMPS 或 GROMACS 的 input 文件就花掉两天改一个参数就得重命名 15 个目录、手动改 42 行 mdp 文件、再逐个检查 log 里有没有Step 0后的报错这不是做科研是在当人肉 Makefile。这个「批量分子动力学模拟脚本和示例」包就是为这种场景生的——它不封装底层引擎不碰 LAMMPS/GROMACS 二进制但用 Python Shell 把建模packmol/tleap、拓扑生成acpype/parmed、任务调度slurm/pbs、日志解析energy/temperature/RMSD 提取、失败自动重试这五层胶水全焊死了。适合正在做构象稳定性筛选、溶剂化自由能粗筛、或需要横向对比力场参数的计算化学实践者。它不是教学 demo是某高校计算平台运维组压箱底的生产级脚本集已支撑过单次超 800 个独立 MD 任务的并发调度。2. 脚本架构与核心模块拆解为什么选 Python Bash 组合而不是纯 Python 或纯 Shell2.1 整体分层设计从输入配置到结果归档的六层流水线整个工具包采用清晰的职责分离输入层YAML 配置文件config.yaml定义体系列表、参数网格、集群资源策略建模层调用 packmol 构建溶剂盒或 tleap 生成 AMBER 拓扑支持 PDB/SMILES 双入口模板层Jinja2 渲染.inLAMMPS或.mdpGROMACS模板变量如{{ temperature }}、{{ nsteps }}全局注入调度层生成 slurm 脚本sbatch.sh并自动sbatch提交支持依赖队列--dependencyafterok:监控层后台轮询jobid状态超时自动 kill 并标记FAILED_TIMEOUT解析层用gmx energy/lammps -log提取关键指标输出 CSVPNG 可视化。提示所有模块通过run_pipeline.py统一入口驱动不强制要求环境变量所有路径均相对project_root/解析避免硬编码路径翻车。2.2 为什么坚持 Python Bash 混合血泪经验告诉你边界在哪纯 Python 调用subprocess.run()管理 HPC 作业看似优雅但实际踩过三个深坑信号传递失效Python 进程被CtrlC中断时子进程如sbatch提交的 job常变成僵尸集群资源卡死日志流阻塞gmx mdrun输出含 ANSI 转义符Python 的stdout.read()会卡在\x1b[?25l隐藏光标处权限继承异常某些 HPC 环境中Python 子进程无法继承module load gromacs/2022.5的环境导致gmx命令未找到。而 Bash 天然适配 HPC 场景trap kill -- -$$ EXIT可靠清理子进程树gmx energy -f edr.xvg 2/dev/null | tail -n 18直接管道解析module load在sbatch.sh中原生生效。我们只用 Python 做高阶逻辑YAML 解析、参数组合、失败分类Bash 做原子操作建模、运行、提取二者通过 JSON 中间文件通信如job_config.json写入当前任务参数Bash 脚本读取后渲染模板。这种切分让每个环节可单独测试——你可以bash scripts/build_box.sh job_001单独验证建模不用启动整条流水线。2.3 核心脚本功能速查表脚本路径语言输入输出关键能力scripts/generate_inputs.pyPythonconfig.yaml,templates/jobs/job_001/{topol.tpr, conf.gro, md.mdp}支持grid_search: {temp: [298,310], conc: [0.1,0.15]}自动生成 4 个任务scripts/submit_job.shBashjob_config.jsonslurm-123456.out,job_001.done自动检测sinfo -p gpuq队列空闲数动态分配--gresgpu:1scripts/parse_log.shBashjob_001/md.logjob_001/summary.csv含step, temp, press, pot_en正则匹配Step.*?Temp.*?Press.*?PotEn跳过WARNING: Listed nonbonded interaction噪声行scripts/plot_summary.pyPythonall_jobs/summary.csvplots/rmsd_vs_temp.png自动按system_name分组用 seaborn 绘制箱线图注意所有 Bash 脚本第一行强制#!/usr/bin/env bash -euxo pipefail-e遇错退出-u禁止未定义变量-x打印执行命令-o pipefail让a \| b \| c中任意环节失败即整体失败——这是 HPC 脚本可靠性的底线。3. 快速上手从零部署到跑通第一个批量任务以 GROMACS 为例3.1 环境准备三步确认你的集群已就绪先别急着 clone用这三条命令验证基础环境是否达标# 1. 确认 GROMACS 可用且版本 ≥ 2021.5因使用 -gpu_id 参数 gmx --version | head -n 1 # 2. 确认 packmol 可调用建模必需 packmol --help | grep Create molecular /dev/null echo ✓ packmol OK || echo ✗ packmol missing # 3. 确认 Slurm 可提交非 PBS/LSF 用户需修改 scripts/submit_job.sh 中的 sbatch 替换逻辑 sinfo -s | grep up /dev/null echo ✓ Slurm OK || echo ✗ Slurm down若任一失败请按错误提示安装对应工具packmol编译需gfortranGROMACS 推荐用 conda 安装gromacs2022.5避免 MPI 版本冲突。3.2 初始化项目结构复制模板并修改最小配置克隆仓库后进入根目录执行# 复制示例配置含 3 个测试体系water_box, protein_ligand, membrane cp -r examples/gromacs_example/* ./ # 修改 config.yaml —— 这是你唯一需要编辑的文件 nano config.yaml打开config.yaml重点修改三处其余保持默认即可# 1. 指定你的体系列表支持 glob 模式 systems: - pdb_files: [examples/pdbs/1abc.pdb] # 蛋白结构 ligand_smiles: c1ccccc1 # 苯将自动用 acpype 生成拓扑 name: benzene_in_1abc # 2. 定义参数扫描网格此处生成 2×24 个任务 parameters: temperature: [298, 310] # K nsteps: 10000 # 步数小规模测试用 dt: 0.002 # ps/step # 3. 集群资源策略根据你的 HPC 调整 cluster: partition: cpuq # 队列名 time: 02:00:00 # 时长 ntasks: 16 # CPU 核数GROMACS 最佳并行度 mem_per_cpu: 4G # 内存逻辑说明systems下每个元素生成一个独立任务目录如jobs/job_001/parameters的笛卡尔积决定该任务下的子任务数如temperature × nsteps生成 4 个md.mdp。name字段会成为输出目录名前缀方便后续grep benzene_in_1abc all_jobs/summary.csv快速筛选。3.3 执行全流程四条命令走完建模→提交→等待→解析# 步骤1生成全部输入文件耗时取决于体系数3 个体系约 15 秒 python scripts/generate_inputs.py # 步骤2提交所有任务自动检测 slurm 队列负载每 5 个任务加 1 秒间隔防 flood bash scripts/submit_all.sh # 步骤3实时监控进度显示 running/pending/done 数量CtrlC 退出 bash scripts/monitor_jobs.sh # 步骤4所有 done 后一键解析全部结果 python scripts/aggregate_results.py执行完aggregate_results.py你会得到all_jobs/summary.csv所有任务的step,temp,pot_en,rmsd时间序列plots/temperature_summary.png按温度分组的 RMSD 箱线图logs/run_pipeline.log完整执行轨迹含每个任务的启动时间、耗时、错误码。参数说明submit_all.sh默认启用--max-concurrent 10同时最多提交 10 个任务避免集群调度器拒绝monitor_jobs.sh每 30 秒轮询一次squeue -u $USER比watch -n 30 squeue更省资源aggregate_results.py会自动跳过job_001/summary.csv缺失的任务只汇总成功项。4. 避坑指南五个真实翻车现场与当场急救方案4.1 现象generate_inputs.py报错ModuleNotFoundError: No module named jinja2原因脚本依赖jinja2渲染模板但未隔离 Python 环境。HPC 登录节点常预装旧版 Python如 3.6而pip install jinja2可能装到用户目录/usr/bin/python3无法加载。解决不要用系统 Python用 conda 创建专用环境conda create -n md_batch python3.9 conda activate md_batch pip install jinja2 PyYAML numpy pandas修改generate_inputs.py第一行#!/usr/bin/env python3→#!/usr/bin/env python3.9或你的 conda 环境路径/path/to/miniconda3/envs/md_batch/bin/python4.2 现象submit_job.sh提交后squeue显示PDpending但 10 分钟不转Rrunning原因Slurm 资源请求超限。例如config.yaml中设ntasks: 128但cpuq队列最大只允许 32 核/任务。解决查看队列限制scontrol show partition cpuq | grep -E (MaxNodes|MaxCPUsPerNode|State)临时降配测试在config.yaml中将ntasks改为16mem_per_cpu改为2G长期方案在scripts/submit_job.sh中添加资源探测逻辑# 获取 cpuq 最大核数 MAX_CORES$(scontrol show partition cpuq | grep MaxCPUsPerNode | awk {print $3}) if [ $NTASKS -gt $MAX_CORES ]; then echo Warning: requested $NTASKS cores queue limit $MAX_CORES, capping... NTASKS$MAX_CORES fi4.3 现象job_001/md.log中出现Fatal error: Too many LINCS warnings任务在 Step 12 就崩溃原因LINCSD 的约束容差constraints h-bonds与初始结构冲突常见于 packmol 生成的溶剂盒未充分能量最小化。解决在templates/em.mdp能量最小化模板中强化约束; 原配置可能只有 500 步改为 2000 步并收紧公差 nsteps 2000 emtol 1000.0 ; kJ/mol/nm原值常为 1000.0此处保持 constraints h-bonds constraint_algorithm lincs lincs_iter 2 ; 增加迭代次数 lincs_order 4 ; 提高精度阶数强制在generate_inputs.py中插入 EM 步骤修改config.yaml的stages字段确保em在md之前执行。4.4 现象parse_log.sh输出的summary.csv中temp列全为0.0原因GROMACS 日志格式变更。新版 GROMACS2022默认关闭tcoupl的详细输出md.log中不再打印Temperature行。解决在templates/md.mdp中显式开启; 在输出控制部分添加 tcoupl v-rescale tc-grps System tau-t 0.1 ref-t 298 ; 关键修复强制输出温度统计 nstcalcenergy 100 ; 每 100 步计算能量 nstenergy 100 ; 每 100 步写入 .edr nstlog 100 ; 每 100 步写入 .log确保 Temperature 行存在或改用gmx energy从.edr提取scripts/parse_log.sh中将grep Temperature md.log替换为gmx energy -f job_001/md.edr -o /tmp/temp.xvg EOF Temperature EOF # 再用 awk 提取最后一列 tail -n 18 /tmp/temp.xvg | awk {print $NF} job_001/temp.csv4.5 现象monitor_jobs.sh显示COMPLETED但job_001/summary.csv为空原因任务虽完成但gmx mdrun因磁盘满No space left on device静默失败.edr未生成parse_log.sh因找不到文件直接退出。解决在submit_job.sh结尾添加磁盘检查# 提交前检查 home 目录剩余空间GROMACS 临时文件巨大 FREE_SPACE$(df ~ | tail -n 1 | awk {print $4}) if [ $FREE_SPACE -lt 10000000 ]; then # 10GB echo ERROR: Less than 10GB free in home, aborting job exit 1 fi在parse_log.sh中增加.edr存在性校验if [ ! -f job_001/md.edr ]; then echo ERROR: md.edr missing for job_001, check mdrun log echo job_001,MISSING_EDR,0,0,0 all_jobs/summary.csv exit 1 fi5. 进阶技巧用失败日志反向定位力场缺陷不止于跑通更要读懂失败5.1 构建「失败模式」数据库把报错文本映射到物理原因批量任务的价值不仅在于成功案例更在于失败样本。我们发现超过 68% 的gmx mdrun崩溃可归为三类物理异常只需解析md.log开头 50 行即可预判日志关键词物理含义应对动作Distance between particles is too small原子重叠建模缺陷回退到packmol阶段增大distance_step参数Potential Energy too large力场参数冲突如电荷过大检查ligand.itp中[ atoms ]的charge值用acpype -c gas重算电荷LINCS WARNING: relative constraint deviation after LINCS键长约束发散积分步长过大降低dt至0.001或改用berendsen压控替代Parrinello-Rahman实现自动化识别只需在scripts/analyze_failures.py中加入import re import pandas as pd def classify_failure(log_path): with open(log_path) as f: head f.readlines()[:50] # 只读前 50 行 log_text .join(head) if re.search(rDistance between particles is too small, log_text): return ATOM_OVERLAP elif re.search(rPotential Energy too large, log_text): return CHARGE_CONFLICT elif re.search(rLINCS WARNING, log_text): return CONSTRAINT_DIVERGE else: return UNKNOWN # 批量分析所有失败任务 failed_jobs [j for j in os.listdir(jobs) if os.path.exists(fjobs/{j}/md.log)] results [] for job in failed_jobs: log_path fjobs/{job}/md.log if Fatal error in open(log_path).read(): results.append({job: job, category: classify_failure(log_path)}) pd.DataFrame(results).to_csv(failure_taxonomy.csv, indexFalse)运行后生成failure_taxonomy.csv你就能一眼看出本次 200 个任务中12 个失败源于原子重叠需优化 packmol 参数3 个源于电荷冲突需重跑 acpype这比盲目重跑全部任务高效十倍。5.2 用「失败任务」反哺建模参数优化一个闭环实验设计某次批量运行中我们发现ligand_smiles: CC(O)O乙酸在 310K 下全部崩溃但 298K 正常。直觉是高温加剧了羧基质子离解而 AMBER99SB-ILDN 力场未覆盖此状态。验证思路从失败任务中提取job_001/conf.gro崩溃前最后结构用pdb2pqr重新质子化pdb2pqr --ff AMBER --ph 7.4 job_001/conf.gro job_001/protonated.pqr用acpype为新结构生成拓扑acpype -i job_001/protonated.pqr -b job_001/protonated将新protonated.top替换原topol.top重新提交job_001。实操细节pdb2pqr需提前conda install -c conda-forge pdb2pqracpype生成的protonated_GMX.top要重命名为topol.top并替换#include路径为避免污染建议在jobs/job_001/下新建revised/子目录存放新文件。5.3 为长期项目建立「参数快照」机制避免复现不可追溯最痛的教训是三个月后想复现某次成功结果却发现config.yaml被改过三次templates/md.mdp也更新了连gromacs --version都升级了。我们的解决方案是每次generate_inputs.py运行时自动生成jobs/job_001/SNAPSHOT/目录包含config_used.yaml当时读取的完整配置template_rendered.mdp实际生成的 mdp含所有 Jinja2 变量值env_info.txtgmx --version,python --version,uname -a在aggregate_results.py中将SNAPSHOT路径写入summary.csv的snapshot_dir列。这样未来任何一行数据都能回溯到精确的软硬件环境。从那以后我每次启动批量任务都强制python scripts/generate_inputs.py --snapshot哪怕多花 2 秒——因为没有快照的 MD 结果就像没有 timestamp 的日志永远无法证伪。希望帮到你。本文还有配套的精品资源点击获取