4. 使用 Python 调用 GROMACS 命令（subprocess）

到了这一步，我们得让 Python 真枪实弹地去指挥 GROMACS 干活了。说白了就是用 Python 的 subprocess 模块启动 GROMACS 的可执行文件，然后拿到结果。这不只是简单的“运行命令”，里面有不少坑。我个人踩过不少，今天跟你掰扯清楚。

4.1 为什么绕不开 subprocess？

你可能会问：“直接在终端里输命令不就行了吗？干嘛要折腾 Python 去调用？”

好问题。想想看，如果你要处理 100 个轨迹文件，每一个都要先跑 gmx trjconv 进行预处理，再跑 gmx rmsf 计算残基波动，最后还要把结果整理成表格。

用手动执行：

• 你得盯着终端，一个一个输命令。
• 输错了要重新来。
• 跑通宵？别想了。
• 出错了也无法自动重试或记录。

用 Python + subprocess 调用：

• 写一个循环，自动遍历所有轨迹文件。
• 遇到错误自动记录日志，跳过或重试。
• 通宵跑完，第二天直接看报表。

我个人习惯把 GROMACS 的调用封装成一个函数，这样在任何脚本里都能直接复用。省心，也少出错。

4.2 subprocess 的基础用法

实现上非常简单，核心就是 subprocess.run()。我贴一段最常用的代码。

import subprocess

# 最基本的调用
result = subprocess.run(
    ['gmx', 'sasa', '-f', 'traj.xtc', '-s', 'topol.tpr', '-o', 'area.xvg'],
    capture_output=True,
    text=True,
    cwd='/path/to/your/workdir'
)

# 检查是否成功
if result.returncode == 0:
    print("SASA 计算完成")
else:
    print("出错了:", result.stderr)

这里要注意几个参数：

• capture_output=True：把标准输出和错误输出都捕获下来，方便我们分析。
• text=True：将输出转为字符串，便于后续处理。
• cwd：指定工作目录。你想想看，如果不指定，Python 会从当前工作目录执行，而你的轨迹文件可能放在别处，路径容易乱。

4.3 捕获并解析输出

很多 GROMACS 工具会在标准输出里打印进度、警告和最终结果。单纯拿到字符串还不够，你得解析出有用信息。

举个例子，gmx rmsf 计算完残基波动后，会打印类似下面的内容：

            Relatively isolated 4.2 Angstrom^2
  Total mean square displacement: 2.3456 nm^2

你需要把这个数字提取出来。我一般用正则表达式来处理，简单粗暴：

import re

def parse_total_msd(output: str) -> float:
    # 匹配 "Total mean square displacement: xxx nm^2"
    pattern = r'Total mean square displacement:\s+([\d.]+)'
    match = re.search(pattern, output)
    if match:
        return float(match.group(1))
    else:
        return None

result = subprocess.run(['gmx', 'rmsf', ...], capture_output=True, text=True)
msd = parse_total_msd(result.stdout)
print(f"总均方根波动: {msd} nm^2")

4.4 错误处理：别让脚本无声无息地崩溃

我刚开始用 subprocess 时，经常忘记检查返回码。脚本跑了一整夜，第二天回来发现中间某个命令报错了，但后面的代码还在傻傻地跑，最后输出的数据全是错的。

基本做法是检查 returncode。但更健壮的方式是主动抛出异常：

def run_gmx(command_list, workdir='.', timeout=3600):
    try:
        result = subprocess.run(
            command_list,
            capture_output=True,
            text=True,
            cwd=workdir,
            timeout=timeout
        )
        result.check_returncode()  # 非零返回码会触发 CalledProcessError
        return result.stdout
    except subprocess.TimeoutExpired:
        raise TimeoutError(f"命令 {command_list[1]} 超时，超过 {timeout} 秒")
    except subprocess.CalledProcessError as e:
        raise RuntimeError(f"GROMACS 命令失败: {e.stderr}")

这样，只要 GROMACS 报错，你的 Python 脚本就会立刻抛异常，并记录下完整的错误信息。你在日志里一看就知道是哪个命令、什么参数、在哪一步出了问题。

4.5 批量执行与并行化

批处理的核心是循环。假设你有 50 个模拟体系，每个体系需要计算 SASA、RMSF 和氢键。我的做法是构建一个任务队列：

如果按顺序跑，50 个体系 × 3 个任务 = 150 个 GROMACS 调用，可能需要几个小时。这时候并行化就很有用了。

Python 里的 concurrent.futures 库是个好帮手。我写个简化示例：

import concurrent.futures
import os

def process_system(sys_id):
    workdir = f'data/{sys_id}'
    # 执行 SASA
    run_gmx(['gmx', 'sasa', '-f', 'traj.xtc', '-s', 'topol.tpr', '-o', 'sasa.xvg'], workdir)
    # 执行 RMSF
    run_gmx(['gmx', 'rmsf', '-f', 'traj.xtc', '-s', 'topol.tpr', '-o', 'rmsf.xvg', '-res'], workdir)
    return sys_id

# 最多同时运行 4 个 GROMACS 实例
with concurrent.futures.ThreadPoolExecutor(max_workers=4) as executor:
    results = list(executor.map(process_system, range(1, 51)))
    print(f"完成 {len(results)} 个体系的处理")

4.6 综合示例：构造一个完整的批处理质控流程

好，我把前面的知识点串起来。下面的脚本会遍历所有子目录，对每个体系执行 SASA、RMSF 和氢键计算，并把关键指标汇总到一个 CSV 文件里。这是我在项目中经常使用的模板。

import subprocess
import re
import csv
import os
from pathlib import Path

def run_gmx(cmd, workdir):
    result = subprocess.run(cmd, capture_output=True, text=True, cwd=workdir)
    result.check_returncode()
    return result.stdout

def parse_sasa(output):
    # 提取总面积，假设格式 "Total area = 123.456 (nm^2)"
    match = re.search(r'Total area\s*=\s*([\d.]+)', output)
    return float(match.group(1)) if match else None

def parse_rmsf(output):
    match = re.search(r'Total mean square displacement:\s*([\d.]+)', output)
    return float(match.group(1)) if match else None

def parse_hbond(output):
    # 提取平均氢键数，假设格式 "Average number of hydrogen bonds: 12.3"
    match = re.search(r'Average number of hydrogen bonds:\s*([\d.]+)', output)
    return float(match.group(1)) if match else None

results = []
base_dir = Path('simulations')
for sim_dir in base_dir.iterdir():
    if not sim_dir.is_dir():
        continue
    sys_id = sim_dir.name
    traj = sim_dir / 'traj.xtc'
    topol = sim_dir / 'topol.tpr'
    if not (traj.exists() and topol.exists()):
        continue

    sasa = run_gmx(['gmx', 'sasa', '-f', str(traj), '-s', str(topol), '-o', 'sasa.xvg'], str(sim_dir))
    rmsf = run_gmx(['gmx', 'rmsf', '-f', str(traj), '-s', str(topol), '-o', 'rmsf.xvg', '-res'], str(sim_dir))
    hbond = run_gmx(['gmx', 'hbond', '-f', str(traj), '-s', str(topol), '-num', 'hbnum.xvg'], str(sim_dir))

    results.append({
        'system': sys_id,
        'sasa': parse_sasa(sasa),
        'rmsf': parse_rmsf(rmsf),
        'hbond': parse_hbond(hbond)
    })

with open('qc_summary.csv', 'w', newline='') as f:
    writer = csv.DictWriter(f, fieldnames=['system', 'sasa', 'rmsf', 'hbond'])
    writer.writeheader()
    writer.writerows(results)
print(f"汇总报告已生成: qc_summary.csv")

这段代码虽然不长，但包含了完整的流程：遍历目录、运行命令、解析结果、输出报告。我建议你把它当作起点，根据实际需求扩展。比如加上错误重试机制，或者记录每次运行的耗时。

好了，这部分内容就这些。你把手里的模拟数据跑一遍，应该就能体会到批处理的威力了。