干这行的人肯定懂一个siesta计算没跑起来多半不是物理模型出了问题而是*.fdf输入文件写得不对。fdf 这东西看着就是个文本但参数、结构、基组、k 点全混在一个文件里手工维护起来非常烦。尤其是结构优化完要更新坐标、批量准备掺杂构型、或者想把参数模板和结构模板分开管理的时候复制粘贴 fdf 段落就成了最磨人的重复劳动。merge-fdf.py就是为解决这个问题写的一个小工具按照 fdf 自身的语法规则把多个 fdf 文件按顺序拼接成一个同名的普通参数和%block结构块以后出现的为准位置保持在第一次出现的地方。脚本只用 Python 标准库不依赖第三方包能跑python3就能用。适合刚接触 siesta、还在手工拼 fdf 的人也适合已经在管理大量模板文件、想把输入文件生成流程化的老手。这篇文章我会把 fdf 的格式拆开讲清楚然后结合merge-fdf.py的完整实现说说脚本背后的设计思路再用三个真实场景演示怎么用。文末还会把我调试过程中踩过的坑和后续扩展方向一并写出来希望对你有帮助。1. 为什么需要 merge-fdf.py三个常见的手拼现场1.1 结构优化后的坐标替换参数没变坐标全变了跑完一轮结构优化下一步通常是把优化后的原子坐标拿回来套用同一套计算参数继续算电子结构、能带或者态密度。这个操作听着简单实际做起来很折腾。我早期是打开输出文件从里面找到最后一步的坐标块然后手动复制替换到输入 fdf 里。坐标少的时候还好几十个原子也能凑合可一旦体系大了或者做了可变晶胞的优化需要同时更新格矢手动操作很容易出错。复制的时候漏一行、多复制一个原子、或者把旧坐标的某几列留下来这类问题几乎没有预兆直到计算跑起来才发现几何构型已经乱了。更麻烦的是坐标块替换完之后还得检查NumberOfAtoms、NumberOfSpecies以及%block ChemicalSpeciesLabel是否同步。原子数不变还好假如你在结构文件里多删了一个原子NumberOfAtoms没同步更新siesta 通常会直接报错或者读入一个错误的结构。这种低级错误很消磨耐心。1.2 批量准备掺杂/吸附构型重复劳动最磨人做表面吸附、掺杂筛选这类工作时经常要准备大量结构只有微小差异的输入文件。比如一个表面模型吸附位点有顶位、桥位、空位三种每种又分不同吸附分子朝向算下来几十个目录。每个目录里的 fdf 参数部分完全一样结构部分只在吸附原子坐标上有区别。如果靠手工复制就要先复制一份模板再改坐标再另存为下一个目录。重复几十次之后人很容易麻木一麻木就会出错。我见过有同事把 A 位点的坐标复制到 B 位点的目录里跑完一个数据点才发现吸附构型复制错了整个计算白跑。这不是细心不细心的问题而是这种流程设计本身就有问题——一旦涉及批量生成输入文件就必须让脚本代替手工去合并、去覆盖、去生成。1.3 模板化参数管理单文件维护最优解还有一种情况不是批量计算而是你想把参数管理得更清晰。一个完整的 fdf 里计算参数、基组定义、k 点设置、结构坐标全都混在一起想单独调整某个环节就得在几百行里找位置。更合理的做法是把 fdf 拆成几个逻辑文件比如params.fdf只管泛函、截断能、收敛判据species.fdf只管ChemicalSpeciesLabel和 PAO 基组structure.fdf只管格矢和坐标。需要调整哪部分就改哪个文件最后用脚本合并成最终输入。这样每个文件都很短可读性高也方便 git 做版本管理。但拆分之后必须有一个可靠的合并工具否则从一个文件变成五个文件反倒给每次提交作业增加了步骤。merge-fdf.py的核心价值就在这里让文件拆分成为一种日常习惯而不是负担。2. 动手之前先拆开 fdf它到底是个什么格式2.1 fdf 的两种基本语法普通键值行和 %block 结构块要写对拼接逻辑第一步是搞清楚 fdf 文件里到底有哪几种语法。Siesta 的 fdf 格式大体上就两类一类是普通键值行一类是%block结构块。普通键值行是最常见的形式一行里包含一个关键词和它的值SystemLabel siesta_run MeshCutoff 400 Ry PAO.EnergyShift 25 meV MD.TypeOfRun CG MD.NumCGsteps 200 MD.MaxForceTol 0.01 eV/Ang这类行遵循关键词 空格 值的规则值里可能带单位。关键词本身大小写不敏感写meshcutoff和写MeshCutoff效果一样但为了可读性通常会保留惯例写法。%block结构块则用于多行数据典型的包括%block LatticeVectors 3.84000 0.00000 0.00000 0.00000 3.84000 0.00000 0.00000 0.00000 3.84000 %endblock LatticeVectors %block AtomicCoordinatesAndAtomicSpecies 0.00000 0.00000 0.00000 1 1.92000 1.92000 1.92000 2 %endblock AtomicCoordinatesAndAtomicSpecies再比如%block ChemicalSpeciesLabel、%block kgrid_Monkhorst_Pack老版本写法、%block PAO.Basis这些都是同一种结构块语法。它们的共同特点是以%block 名字开头以%endblock结尾中间是多行内容。这里必须注意%endblock后面一般可以跟名字也可以不跟但写全名字会让文件更清晰也方便脚本做校验。fdf 格式的核心语法我用一张表总结一下语法类型示例用途普通键值行MeshCutoff 400 Ry单个标量参数%block 结构块%block LatticeVectors ... %endblock格矢、坐标、基组等多行数据注释行# 这是注释说明性文本会被解析器忽略空行无分隔段落不影响语义2.2 拼接最容易出错的两个地方第一是 block 的多行属性。很多人第一次写合并脚本时会天真地按行去拼接把两个文件的每一行按顺序堆到一起。这在小样本下可能碰巧能跑通但一旦遇到%block结构块按行拼接就会把块内容拆得七零八落甚至把一个文件的后半个 block 和另一个文件的 block 残留拼在一起。真正稳妥的做法是先把文件切成块级别的内容单元再对单元做合并。第二是结构定义的上下文关系。坐标块%block AtomicCoordinatesAndAtomicSpecies里的原子序号必须和%block ChemicalSpeciesLabel里定义的物种序号一一对应。如果你合并时只更新了坐标块却忘了同步物种定义计算会直接报错。同样如果你的格矢可变优化后的LatticeVectors也要一并替换进去否则会出现坐标和格矢不匹配的物理问题。2.3 一个关键认知fdf 的后者覆盖前者规则fdf 文件在解析的时候对于同一个关键词多次出现的处理逻辑一般是后面出现的覆盖前面出现的。这一点很关键因为它直接影响合并脚本的设计方向。如果你的某个文件里可能出现同一个参数的多个定义最安全的合并策略就是顺着文件顺序逐个处理后遇到的同名参数替换之前的同名参数。这样既能保证无脑拼接时可能产生的重复定义问题被消除又能很好地适配用后面的结构文件覆盖前面的模板文件这种典型用法。因此merge-fdf.py的核心语义可以定成合并后同一个普通键或同一个 block 只保留一个以所有输入文件中最后一个出现的为准。3. merge-fdf.py 的核心实现先把文件切成块再按覆盖规则组装3.1 整体设计思路解析、合并、写出三步走脚本整体分为三个层次解析层、合并层、输出层。解析层负责把单个 fdf 文件解析成一个有序的条目列表。每个条目可能是普通键值行、%block结构块、或者注释/空行。解析层不关心数值是什么只负责把结构完整地切分出来。合并层遍历所有输入文件维护一个当前合并结果列表。遇到普通键或 block 时先看是否已经出现过同名条目如果出现过就在原位置替换内容如果没出现过就追加到尾部。注释和空行保持原样追加保证文件读起来不至于太奇怪。输出层把合并结果写回文本文件统一使用 Unix 换行符避免 Windows 和 Linux 混用时的格式问题。3.2 fdf 解析器识别 block 和普通键解析器的核心是一个顺序扫描循环用正则判断当前行是%block开始、%endblock结束、还是普通键值行。这里我把普通键值的单位部分也一并保留在原始行里不拆分这样后续写出去的时候不会因为重新拼接而丢失单位。import re import sys def parse_fdf_file(filepath): 把单个 fdf 文件解析成条目列表。 每个条目是一个元组 (kind, name, content, lineno) kind 取值 block - %block 结构块 key - 普通键值行 raw - 注释或空行 content 保留原始文本内容方便原样输出。 entries [] with open(filepath, r, encodingutf-8, errorsreplace) as fh: lines fh.readlines() i 0 n len(lines) while i n: line lines[i] stripped line.strip() # 空行和注释作为 raw 条目保留 if not stripped or stripped.startswith(#): entries.append((raw, , line, i 1)) i 1 continue # 匹配 %block 开头 block_match re.match(r^\s*%block\s(\S)\s*$, line, re.IGNORECASE) if block_match: block_name block_match.group(1).lower() buf [line] i 1 closed False while i n: buf.append(lines[i]) if re.match(r^\s*%endblock, lines[i], re.IGNORECASE): closed True i 1 break i 1 if not closed: print(f警告: {filepath} 中的 %block {block_name} 没有闭合, filesys.stderr) entries.append((block, block_name, .join(buf), block)) continue # 普通键值行第一个空白字符之前的字段当作 key parts stripped.split(None, 1) if len(parts) 1: key parts[0].lower() entries.append((key, key, line, i 1)) else: entries.append((raw, , line, i 1)) i 1 return entries注意几个细节。第一所有 key 和 block 名都统一转成了小写。因为 fdf 关键词本身大小写不敏感用小写做指纹可以避免LatticeVectors和latticevectors被当成两个不同条目。第二遇到%block时用循环把整个块收集起来而不是只记一行。这样后续替换时才能保持块内容的完整性。第三对于未闭合的 block脚本会打出一条警告。这个分支看起来不起眼但实际很有用。我调试的时候遇到过 fdf 文件被截断的情况如果没有这个发现合并出来的文件会让 siesta 在解析阶段直接崩溃排查起来很痛苦。3.3 覆盖合并逻辑位置不变内容替换解析之后合并的核心就是一个有记忆的顺序遍历。我用一个列表保存最终条目用一个字典记录每个 key/block 在列表中的位置。这样遍历到同名的新条目时可以直接替换掉旧条目同时保持旧条目的原始位置。为什么保持原始位置很重要因为 fdf 文件里有些参数之间存在上下文习惯。比如AtomicCoordinatesFormat Ang通常写在坐标块附近如果你把所有结构块都挪到文件末尾虽然大部分情况下 siesta 也能读但读文件的人会觉得很奇怪。保持位置不变实际效果就是原位更新和手动编辑文件的直觉一致。def merge_files(input_files): 按顺序合并多个 fdf 文件返回合并后的条目列表。 merged [] index {} for filepath in input_files: for ent in parse_fdf_file(filepath): kind, name, content, _ ent[0], ent[1], ent[2], ent[3] if kind raw: # 注释、空行直接保留 merged.append(ent) continue # block 和 key 用不同的前缀做指纹防止互相覆盖 fingerprint (b: if kind block else k:) name if fingerprint in index: pos index[fingerprint] merged[pos] ent # 替换旧条目位置不变 else: index[fingerprint] len(merged) merged.append(ent) return merged def write_fdf(entries, out_path): 把条目列表原样写回文件统一使用 Unix 换行。 with open(out_path, w, encodingutf-8, newline\n) as fh: for ent in entries: fh.write(ent[2])这十几行就是合并逻辑的全部。它没有做任何数值层面上的处理因为拼接 fdf 的底线是不改变原文的含义只做结构上的重组。3.4 命令行入口一个文件搞定不依赖第三方库为了让脚本在命令行下用起来顺手我再加一个简单的参数入口。支持指定多个输入文件和一个输出文件并增加一个--check选项来做基础校验。import argparse def check_consistency(entries): 检查合并结果中 NumberOfAtoms 与坐标块行数是否一致。 natoms None coord_lines None for kind, name, content, _ in entries: if kind key and name numberofatoms: parts content.split(None, 1) try: natoms int(parts[1]) except (IndexError, ValueError): natoms None if kind block and name atomiccoordinatesandatomicspecies: lines [ln for ln in content.splitlines() if ln.strip() and not ln.strip().startswith(%)] coord_lines len(lines) if natoms is not None and coord_lines is not None and natoms ! coord_lines: print(f警告: NumberOfAtoms 为 {natoms}坐标块却有 {coord_lines} 行, filesys.stderr) return False return True def main(): parser argparse.ArgumentParser( descriptionmerge-fdf.py: 合并多个 Siesta fdf 文件 ) parser.add_argument(inputs, nargs, help输入 fdf 文件按从左到右的顺序合并) parser.add_argument(-o, --output, defaultmerged.fdf, help输出文件路径默认 merged.fdf) parser.add_argument(--check, actionstore_true, help合并后检查原子数与坐标行数是否一致) args parser.parse_args() entries merge_files(args.inputs) if args.check: check_consistency(entries) write_fdf(entries, args.output) print(f已生成 {args.output}) if __name__ __main__: main()至此merge-fdf.py从解析到输出就完整了。整份代码不到一百行没有什么高深算法但它在实际工作流里的用处非常大。4. 实战演示三种场景下 merge-fdf.py 的具体用法4.1 场景一结构优化后更新坐标假设你有一个base.fdf里面是完整的计算参数和初始结构。优化结束后你想保留所有参数只把坐标和格矢换成优化后的最终结构。这时你可以先把最后一帧结构单独抽成一个final_struct.fdf。比较省事的方式是用 ase 读取 Siesta 的轨迹文件把最后一帧原子坐标写出来再拼成一个干净的结构块。下面是示例代码from ase.io import read atoms read(siesta.MD, index-1) with open(final_struct.fdf, w) as f: f.write(AtomicCoordinatesFormat Ang\n) f.write(%block AtomicCoordinatesAndAtomicSpecies\n) for atom in atoms: # 注意: species 编号必须与 ChemicalSpeciesLabel 一致 sp atom.tag if atom.tag else 1 pos atom.position f.write(%20.10f %20.10f %20.10f %d\n % (pos[0], pos[1], pos[2], sp)) f.write(%endblock AtomicCoordinatesAndAtomicSpecies\n)如果你的优化过程允许晶胞变化还需要把最后的LatticeVectors也提取出来。很多版本的siesta.MD里会周期性打印格矢提取时注意取最后一帧对应的格矢不要取初始格矢。然后执行python3 merge-fdf.py -o run.fdf base.fdf final_struct.fdf这样生成的run.fdf参数部分完全继承自base.fdf坐标部分自动替换成最后一帧结构所有结构块的原始位置保持不变。执行完之后建议加一个--checkpython3 merge-fdf.py -o run.fdf base.fdf final_struct.fdf --check这个检查会统计NumberOfAtoms和坐标块行数万一你在提取结构时漏了一个原子它会在提交作业之前就提醒你。4.2 场景二批量生成吸附构型输入文件批量场景下merge 的优势更加明显。你只需要维护好两个文件夹一个是存放公共参数模板的目录一个是存放各构型结构文件的目录。比如结构文件命名为site1.fdf、site2.fdf、site3.fdf每个文件里只写坐标块和格矢块for site in site1 site2 site3; do python3 merge-fdf.py -o run_${site}.fdf params.fdf ${site}.fdf done循环跑完run_site1.fdf、run_site2.fdf、run_site3.fdf就都生成好了。每个文件里参数部分完全一样结构部分来自对应的构型文件。接下来不管是本地批量跑还是丢到集群都可以接一个作业生成脚本继续往下走。这里有一个需要注意的地方如果某个构型文件里没有写NumberOfAtoms那么合并后的文件会沿用params.fdf里的NumberOfAtoms。如果你的几个构型原子数不同就必须在各自的结构文件里显式写入NumberOfAtoms否则合并结果就是错的。这个容易出现我建议在每个结构文件头部都写上原子数和物种数不要依赖模板去猜。4.3 场景三参数模板与结构模板分离管理我现在的日常工作方式是给每个项目建一个templates/目录里面放几类模板文件作用params.fdf泛函、截断能、收敛判据、MD 参数等计算设置species.fdfChemicalSpeciesLabel与 PAO 基组定义kpoints.fdfk 点采样设置cell.fdf格矢和原子坐标结构定义需要生成一个新作业时按需合并python3 merge-fdf.py -o newjob.fdf params.fdf species.fdf kpoints.fdf cell.fdf这样做的收益是参数调整只动params.fdf一个文件结构替换只动cell.fdf一个文件其他人接手时也能一眼看清整个计算设置的分层结构。配合 git 之后每次计算的输入文件从哪来、改过什么都清清楚楚。5. 实测踩坑从能跑到跑对5.1 编码和换行符Windows 上打开 Linux 文件容易翻车siesta 主要跑在 Linux 集群上但很多时候模板文件是在 Windows 上编辑的。Windows 的文本编辑器和 Linux 的换行符不一致文件里会混入\r\n。如果不做处理merge 生成的 fdf 可能在某些环境下解析出奇怪问题。我在write_fdf里显式指定了newline\n读取时用了encodingutf-8, errorsreplace。读取时用errorsreplace可以避免遇到非法编码时直接崩溃。这个小细节带来的好处是无论输入文件是 Windows 还是 Linux 格式输出文件统一是 Unix 换行提交到集群上不会因为换行符问题报错。需要说明的是我并没有在读取时把\r\n显式转成\n因为 Python 打开文本文件时会做通用的换行转换只要write_fdf输出时统一成\n就够了。实际用下来这个方案最省心。5.2 坐标精度不能顺手丢写解析器时一个很容易犯的错误是对坐标数值做格式化。比如读取坐标块里的每一行把字符串转成 float 再拼回去。浮点数在打印时如果位数不够会悄悄丢失精度。对原子位置来说哪怕丢掉一位有效数字优化可能就白做了。所以我一开始就没打算在解析阶段去理解数值的内容。整个脚本把普通键值和 block 的内容都当作透明文本处理原样保留、原样输出。坐标位数在用户自己的结构文件里是多少合并出来就是多少脚本绝不碰数值。这是一个很重要的设计取舍拼接脚本只负责结构和去重不负责格式化数据。5.3 同名键在不同位置覆盖的陷阱后出现的同名条目会替换先出现的同名条目这个规则大多数时候很好用但也会带来一个隐蔽的问题。比如params.fdf里你已经写好了SystemLabel run_001后面的kpoints.fdf里不小心也写了一个SystemLabel run_kpt合并之后输出文件里的系统标签就变成了run_kpt作业输出前缀全变了。这种问题不算致命但会浪费一轮排查时间。我的对策是两层第一所有模板文件保持简洁最上层的公共参数只在一个文件里出现第二在用脚本生成最终输入后打开文件扫一眼开头和结尾确认没有明显不符合预期的重复定义。如果文件很多也可以用grep -i systemlabel快速确认。5.4 合并后必做的两道检查第一道检查是差异对比。可以把合并结果和原始的 base 文件做一次 diff看看到底哪些行被替换了、哪些行被追加了。特别是在更新结构坐标时diff 能很明显地暴露坐标块是否被正确替换以及NumberOfAtoms是否出现了重复定义。第二道检查就是前面提到的--check。这个检查不复杂但对吸附构型这类原子数可能变化的场景非常管用。我建议在任何更新过结构的合并操作中都加上--check成本几乎为零却能提前发现最典型的低级错误。6. 扩展思路让它长成你自己的小工具6.1 给脚本加一个 --diff 模式实际用起来之后你会发现合并前看差异和合并后看差异是两种完全不同的需求。有些模板文件很长你只想知道两个版本之间到底改了什么。可以在脚本里增加一个--diff参数直接调用系统diff命令对比两个合并结果。更简单的方式其实是把merge-fdf.py和命令行 diff 组合使用python3 merge-fdf.py -o old.fdf params_old.fdf cell.fdf python3 merge-fdf.py -o new.fdf params_new.fdf cell.fdf diff old.fdf new.fdf这样能快速看到参数变更对整体输入的影响适合在做收敛性测试时使用。6.2 反向需求从一个 fdf 里单独抽出结构块拼接有了反向提取也是常用功能。有时候你从别人的项目里拿到一个完整的 fdf想把里面的结构部分抽出来作为自己的模板配合参数文件重新组织项目。这个需求不需要再写一个独立工具只要用parse_fdf_file把 block 挑出来from merge_fdf import parse_fdf_file for kind, name, content, _ in parse_fdf_file(example.fdf): if kind block and name in (latticevectors, atomiccoordinatesandatomicspecies): print(content)有了解析层这类小需求都可以顺手实现不用再从零写正则。6.3 与自动化作业生成流水线结合最后说一个我个人的工作流习惯。我现在不会手动运行merge-fdf.py而是把它嵌到整个作业生成流水线里。大体流程是先用 Python 脚本批量生成所有结构文件再用循环调用merge-fdf.py把每个结构的 fdf 拼好接着生成对应的 PBS/SLURM 作业脚本最后统一提交。这样的好处是每个目录下的输入文件都可以随时从模板重新生成。即使某个目录被误删了只要跑一遍生成脚本所有输入就能恢复原样。对我来说merge-fdf.py不再是一个孤立的脚本而是整个计算流程中一个可靠的基础设施。如果你也经常和 fdf 文件打交道建议从最基础的合并开始用慢慢加上适合自己项目的参数和扩展。这个脚本不大但它能把那些重复、容易出错的手工操作从日常工作中彻底去掉省下来的时间足够你多跑好几个构型了。