📚 课程目录第1章计算化学任务自动化工程导论

1. 计算化学任务自动化工程导论

计算化学的日常工作，往往不是“想一个理论”，而是“把同一套理论在几十上百个分子上重复执行”。一次几何优化、一次频率计算、一次溶剂化单点能，单独看都不复杂；但当体系从 1 个变成 200 个，从一种泛函变成四种泛函，从气相变成三种溶剂模型时，真正消耗精力的就不再是化学判断，而是文件命名、输入生成、任务排队、结果解析和失败重跑。自动化工程要解决的，正是这类“重复但必须可靠”的问题。

可以把自动化在计算化学中的角色理解为三层递进。第一层是脚本化：用一段 Shell 或 Python 把手工敲的命令串起来，例如批量生成 Gaussian 输入文件、循环调用 g16、把 .log 里的能量抓出来。脚本化的核心价值是省去重复输入，但它通常假设“一切顺利”，一旦某个任务因 SCF 不收敛或磁盘写满而失败，脚本可能继续往下跑，留下难以追溯的残缺结果。第二层是工作流化：把任务拆成有明确依赖关系的步骤，例如“结构准备 → 几何优化 → 频率验证 → 单点能 → 结果汇总”，并让每一步的输入来自上一步的产物、失败时能停止或重试。工作流化关注的是步骤之间的数据契约与执行顺序。第三层是工程化：在工作流之上加入版本管理、参数记录、日志、校验、可复现的环境与失败恢复机制，使整套流程不仅能跑通一次，还能被别人复现、被自己半年后重跑、在集群上稳定扩展。

三层不是互斥的替代关系，而是能力叠加：脚本化解决“少打字”，工作流化解决“按顺序、有依赖”，工程化解决“可复现、可维护、可扩展”。判断一个计算化学项目是否需要工程化，关键看它是否要重复执行、是否要交付给他人、是否要在参数变化下比较结果。

要理解自动化收益从何而来，需要先看清一个典型计算化学任务的生命周期。以“对一组候选分子做气相几何优化并比较相对能量”为例，其生命周期大致包含：定义体系与计算级别（泛函、基组、溶剂、电荷与自旋多重度）、生成初始结构、构造输入文件、提交与调度、监控运行状态、解析输出、校验收敛与频率、汇总数据、归档与记录。这个链条中，真正需要化学家判断的是“计算级别是否合理”“结果是否物理可信”；而生成输入、提交、解析、汇总这些环节，恰恰是最适合自动化的部分。

上图中，蓝色主链是任务从定义到运行的自然顺序，绿色支链是结果处理与归档，橙色虚线框标出自动化可以介入的环节。可以看到，自动化并不是替代化学判断，而是把“定义级别”之后的重复劳动接管下来，并把“校验与判据”变成可编程的检查点。

下面给出一个最小但完整的工程化示例，用 Python 组织一次批量几何优化。它不依赖任何特定计算软件的高级接口，只做三件事：从模板生成输入、按命名规则提交、解析输出中的能量与收敛状态。示例参数为教学假设，实际使用时需替换为本地软件路径与真实分子。

# batch_opt.py —— 教学示例，参数为假设值
import os, subprocess, re, csv

MOLECULES = ["mol_01.xyz", "mol_02.xyz", "mol_03.xyz"]
LEVEL = {"functional": "B3LYP", "basis": "6-31G(d)", "charge": 0, "mult": 1}
TEMPLATE = """%chk={name}.chk
#p {functional}/{basis} opt freq

{name}

{charge} {mult}
{coords}
"""

def build_input(xyz_path):
    name = os.path.splitext(os.path.basename(xyz_path))[0]
    with open(xyz_path) as f:
        lines = f.read().strip().splitlines()
    coords = "\n".join(lines[2:])          # 跳过原子数与注释行
    return name, TEMPLATE.format(name=name, coords=coords, **LEVEL)

def parse_log(log_path):
    text = open(log_path, errors="ignore").read()
    energy = re.findall(r"SCF Done:\s+E\([^)]+\)\s+=\s+(-?\d+\.\d+)", text)
    converged = "Stationary point found" in text
    freqs = re.findall(r"Frequencies --\s+(-?\d+\.\d+)", text)
    return (float(energy[-1]) if energy else None, converged,
            min(float(f) for f in freqs) if freqs else None)

rows = []
for xyz in MOLECULES:
    name, inp = build_input(xyz)
    with open(f"{name}.gjf", "w") as f:
        f.write(inp)
    subprocess.run(["g16", f"{name}.gjf"], check=False)   # 实际需按集群调度改写
    e, conv, fmin = parse_log(f"{name}.log")
    rows.append({"name": name, "energy_Hartree": e,
                 "converged": conv, "min_freq_cm-1": fmin})

with open("summary.csv", "w", newline="") as f:
    writer = csv.DictWriter(f, fieldnames=rows[0].keys())
    writer.writeheader(); writer.writerows(rows)

这段代码已经具备脚本化的骨架，但离工程化还有距离。它缺少失败重试、缺少对 g16 返回码的判断、缺少参数与软件版本的记录，也没有把“虚频为负”作为失败条件显式处理。一个更工程化的版本至少应做到：每个任务有独立目录与日志；提交前检查输入文件是否完整；运行后先判断进程返回码，再判断 SCF 是否收敛、频率是否全为正；对失败任务记录原因并支持单独重跑；最后把计算级别、软件版本、提交时间写入汇总表。这样，当三个月后需要补充一个分子时，你不需要回忆当时用了什么基组。

判断自动化是否“够用”的一个实用标准：如果明天换一台机器、换一个合作者，他能否仅凭你留下的脚本、参数记录和目录结构，复现出同样的相对能量排序？能，就达到了工程化的基本要求。

需要强调的是，自动化工程并不追求一开始就搭建庞大框架。更合理的路径是：先用脚本消除最痛的重复劳动，再把脚本中反复出现的依赖关系抽成工作流，最后针对“必须复现、必须交付”的部分补上记录与校验。对于一次性探索，脚本化足矣；对于要写入论文支撑材料、要被审稿人追问参数的计算，工程化就是必要投入。计算化学的自动化，本质上是用工程手段保护化学判断的可靠性——让机器负责重复与记录，让人负责选择与解释。