简介：面向材料计算与凝聚态物理研究者的 DFT 预处理与后处理开源工具箱，以 Python 模块形式封装，旨在简化 QE、Abinit、ELK 等主流 DFT 代码的输入文件构建和输出结果分析，相较于直接阅读官方文档或逐个调整变量，这一方式能帮助新用户更快上手。设计理念是用户只需提供少量参数，即可快速生成计算输入，并对能带、态密度、投影态密度等结果进行读取、分析与可视化，显著降低初学者面对大量变量的学习成本。压缩包共 329 个文件，约 21.99MB，核心为 30 个 Python 脚本，同时包含 png 示例图、in/kpf 输入文件、out/dat/npz 结果文件、paw/upf 赝势文件以及 abinit/QE 相关的 fatbands 与 pdos 数据，目录结构完整，适合作为 DFT 流程学习和二次开发的参考。已有 430 人浏览学习。通过该工具箱可快速搭建典型计算任务，理解不同代码的输入输出差异，并直接复用或修改示例脚本用于自己的科研项目。

1. DFTtoolbox：MATLAB 里搭 DFT 前处理与后处理的脚手架

做材料计算的人总有一段时间被输入文件折磨：换个结构就要重排坐标，改个精度就要翻文档找关键字，跑完又得用脚本去抠能带、态密度和总能量。DFTtoolbox 是一个基于 MATLAB 的开源工具箱，专门干两件事：把结构参数快速打包成主流 DFT 程序能识别的输入文件，再把跑完的输出文件解析成语义清晰的 MATLAB 变量。这套源代码适合核心问题是物理结果、不打算为格式细节反复折腾的从业者，也适合带学生做课题时减少低级格式错误。之前在某公司做材料模拟时，我几乎每周都要处理几组新结构，这个工具箱让我把精力留在看结果而不是改文件。

2. 输入文件生成：从结构到计算卡，让 DFT 前处理告别手写硬编码

2.1 手工输入卡的痛点与工具箱的定位

DFT 计算的第一步是把晶格常数、原子坐标、截断能、k 点网格、交换关联泛函这些信息写进一个程序能解析的文本文件。以某主流平面波程序为例，至少需要四类文件：结构文件描述晶格和原子位置，参数文件控制计算类型与精度，k 点文件决定布里渊区采样密度，另外还要指定赝势文件和原子种类。手工去改这些文件有个老生常谈的问题——结构一换，坐标就得重排，分数坐标和笛卡尔坐标混用时更是重灾区。

常见做法是给每个结构写一份模板，用脚本填充变量。但模板散落在不同项目里，很容易出现过期关键字、单位不统一、赝势路径写死之类的问题。DFTtoolbox 的定位就是把这些模板收拢成一个配置化入口：你在 MATLAB 里用一个结构体描述体系，再通过函数设置计算参数，最后调一个写入函数，所有文件就按约定格式落到目录里。

选择 MATLAB 而不是纯文本脚本，主要是许多人做材料数据处理的主工作台就是 MATLAB，结构分析、可视化、后续的机器学习扩展都能延续同一套环境。而且 MATLAB 的 struct 和 table 在处理异构数据上比正则表达式直观得多；工具箱把读文件、解析关键字、构建输出这些重复劳动封装成有返回值的函数，调试时可以直接看变量内容，不用盯着抽日志的正则发呆。当然 Python 生态也很成熟，但如果你整个分析流程已经在 MATLAB 里，多引入一个工具箱比重新搭一套脚本更划算。

2.2 结构体、赝势目录与 k 点网格的默认约定

我一般会把晶格信息、原子种类、坐标和坐标类型都放在一个结构体里。工具箱里常见的约定是：a、b、c 是晶格长度，alpha、beta、gamma 是角度，cell_type 决定晶系的归类，coordinates 按分数坐标给出，species 用字符串数组顺次对应每个原子。坐标的约定直接影响对称性识别，所以我会默认全部用分数坐标，只有在做分子或者孤立体系时才切成笛卡尔坐标。

赝势目录在工具箱里是一个配置项，默认指向一个集中存放赝势的路径，每个元素种类按元素符号查找。集中管理的好处是换机器、换结构时不用逐个修改路径。k 点网格的默认逻辑是：金属体系推荐 Gamma-centered 网格，绝缘体可以用 Monkhorst-Pack；对于六方格子，k 点在 a、b 方向取 9×9，c 方向按倒空间比例适当压缩到 3 或 5。不是每个体系都适合同一套默认值，但有个可复用的起点总比每次翻手册强。

参数文件里最容易出错的不是物理量本身，而是单位。平面波程序里截断能与能量窗口的单位通常是电子伏特，但有些程序会接受其他单位制，差一个换算系数就是白跑一趟。工具箱的输入函数在设置项后面都会校验单位标志，这一层静态检查帮我拦住过不少低级错误。

2.3 一段最小可用的输入生成流程

下面这段代码演示了从结构体到输入文件的最小流程：

% 构建晶格与原子结构，六方晶系，三层原子 cryst = struct('a', 3.85, 'b', 3.85, 'c', 9.42, ... 'alpha', 90, 'beta', 90, 'gamma', 120, ... 'cell_type', 'hexagonal'); cryst.species = {'Ti', 'O', 'O'}; cryst.coords = [0.3333, 0.6667, 0.25; ... 0.6667, 0.3333, 0.75; ... 0.5, 0.5, 0.5]; % 生成输入文件，写到 calc01 目录 dft_input = dfttoolbox.input.vasp(cryst, ... 'encut', 520, 'kpoints', [9 9 3], ... 'xc', 'PBE', 'sigma', 0.1, 'ismear', 1); dft_input.write('calc01');

这段代码的流程是：先定义晶格和原子坐标，再调用 vasp 输入构造函数，返回一个输入对象，最后调用 write 方法把结构化文件写入目录。整个过程没有碰任何文本模板，格式细节由工具箱统一处理。

参数的物理含义需要交代一下：encut 是平面波截断能，单位 eV，决定基组完备程度；kpoints 是 [kx, ky, kz] 的 Gamma 中心网格密度；xc 指定交换关联泛函，PBE 是 GGA 类的常用选择；sigma 和 ismear 是展宽方式与展宽宽度，金属体系推荐 ismear=1、sigma=0.1 eV，绝缘体建议改用 ismear=-5 并忽略 sigma。第一次跑某个新材料时，我通常先按上面这套参数走一遍，再根据结果调精度，而不是直接上高精度硬算。

传入函数的关键参数汇总如下，方便你根据自己的体系改。

公共参数确认后，不同结构只需要变更 cryst 内容，其他逻辑完全复用。这也是我后来习惯用这个工具箱而不是保留一堆模板脚本的原因——结构变了，代码不用改。

3. 结果解析：能带、态密度与能量曲线的 MATLAB 读法

3.1 能带数据读取与高对称路径对齐

计算跑完后，能带数据通常分散在输出文件里，一个常见坑是不同 DFT 程序能带文件格式差异很大。以某主流程序的输出为例，能带文件的头部是元数据，后面每一块对应一条 k 点路径，块内每行是一个本征态的能量。手工解析的问题在于：头部行数、态数量标记、以及是否包含费米能级在不同版本里不一致。

工具箱的做法是先扫描文件头确定 nbands、nspin 和 k 点数，再按块读取，同时把费米能级从日志文件里解析出来。读取之后返回的 bands 对象包含 eigs、kpath 和 fermi 三个字段，eigs 是 nk×nband 的矩阵，kpath 是逐条 k 点对应的倒空间距离。画能带图时，横轴需要把 k 点坐标按高对称点距离换算成累计长度，否则形状会失真。

% 读取能带数据，自动定位费米能级 bands = dfttoolbox.result.read_bands('calc01/EIGENVAL', ... 'fermi', 'auto', 'nbands', 96, 'nspin', 1); % 沿高对称路径画出能带，能量零点取费米能级 figure; dfttoolbox.plot.bandstructure(bands, ... 'highsym', {'Gamma', 'M', 'K', 'Gamma'}, ... 'ylim', [-4, 4], 'color', [0 0.45 0.74], ... 'linewidth', 1.6);

代码逻辑：第一行读取能带并自动寻找费米能级；第二行画图时指定高对称点序列和能量窗口。fermi 设为 auto 时，工具箱会从日志里自动解析费米能级，不必手动填，前提是日志文件没有被截断。

高对称路径并非随意指定。六方晶系常用 Gamma-M-K-Gamma 循环，立方晶系常用 Gamma-X-W-K。工具箱内置了几种常见格子的对称点坐标，画图时只需提供对称点名序列，它会自动匹配倒空间坐标并换算累计长度。

注意 nbands 需要包含全部价带和若干导带。一般用电子总数除以 2 再预留一部分空带，预留量太小会导致高能级处看到截断而非真实物理。我曾经图省事把 nbands 设成刚好等于价带数，结果导带底部图像全被切断，重新跑了一遍才确认不是物理问题。体系结构、原子数目变化时，电子数会变，对应的 nbands 也要跟着更新，最稳妥的做法是直接读取输出头部给出的最大带数，而不是硬编码。

3.2 态密度重采样、费米能级平移与分波投影

态密度部分通常要处理展宽和能量网格。输出文件里会给出一个能量网格，但不同段的网格间隔可能不一致，尤其是能量窗口比较宽时，直接绘图会出现疏密不匀。工具箱读取态密度后，会按设定的能量区间做线性重采样，确保每个点间隔相等，方便后续积分和对比。

费米能级平移也是关键一步。态密度一律以费米能级为零点展示，才能直观判断哪个峰在费米面附近。我处理 DOS 数据时习惯先读取原始数据，再统一减去费米能级，之后所有图表都用同一个约定。

% 读取态密度，能量范围为 -20 到 20 eV，重采样 801 点 dos = dfttoolbox.result.read_dos('calc01/DOSCAR', ... 'emin', -20, 'emax', 20, 'nedos', 801); % 费米能级平移到零点 dos.e = dos.e - dos.fermi; % 画出总态密度和两条分波态密度 dfttoolbox.plot.dos(dos, 'components', {'Ti-3d', 'O-2p'}, ... 'xlim', [-8, 6], 'legend', true);

读出来的结构里，projection 字段会随计算配置变化。若计算时开启了投影选项，才会有分波数据，画分波 DOS 前要确认该字段存在。另一个容易被忽略的问题是展宽参数：读取态密度时工具箱默认用高斯展宽做重采样，如果原输出文件已经用展宽处理过，重采样后可能把峰压低。这种情况可以把展宽参数调小，或者直接用原始能量网格。

常见误用是把总态密度和分波态密度混在一起比较数值大小。总 DOS 的单位是 states/eV，而分波 DOS 是投影到单个轨道上的贡献，两者数值量级不同，直接堆叠会让分波曲线几乎看不见。画的时候给分波曲线单独用右侧坐标轴，或者先归一化再叠加。

3.3 总能量与磁矩收敛曲线的提取

能量收敛判断是每个计算都要做的检查。输出文件里，自由能那行是每个离子步或电子步的总能量，混合了不同收敛层级的数据，直接提取会得到一条抖动的曲线。工具箱读取时按标记行分别归类，把电子步和离子步分开，这样就不会把中间电子步的能量当成未收敛的信号。

% 从输出中读取分步能量与磁矩 conv = dfttoolbox.result.read_convergence('calc01/OUTCAR'); % 检查最后两步离子步的能量差，单位换算成 meV/atom n_atoms = length(cryst.species); dE_last = diff(conv.energy_ion(end-1:end)) / n_atoms * 1000; fprintf('最后两步离子步能量差: %.2f meV/atom\n', dE_last); % 绘制能量-离子步曲线，观察收敛趋势 figure; plot(conv.step_ion, conv.energy_ion, 'o-'); xlabel('离子步'); ylabel('总能量 (eV)'); title('总能量收敛曲线');

这里的 conv 对象包含 energy_ion、step_ion 等字段。energy_ion 是每个离子步结束后输出的自由能，step_ion 是对应的步号。能量差小于 1 meV/atom 时我才会继续后续属性计算，否则就回去检查截断能和 k 点网格。

磁矩曲线对磁性材料特别有用。如果输出里每个离子步附带了磁矩信息，工具箱会把磁矩和能量放在同一个结构里，方便观察结构弛豫过程中磁矩有没有跳变。磁矩在某个步数突然翻转，通常说明初始磁矩设置离稳态太远，需要调整初始磁矩再重新跑，否则后面能带和态密度的自洽结果都不可信。

读这些结果时，尽量用脚本把数据全部抓成 MATLAB 变量并保存到 .mat 文件，二次分析就从 .mat 里加载，避免反复读取大文本。输出文件在 MATLAB 里逐行扫描很慢，这一层缓存能显著提速。

4. 参数怎么调：截断能、k 点密度与对称性的三角平衡

4.1 截断能与 k 点密度的搭配原则

DFT 结果有用的前提是数值收敛。截断能决定平面波基函数数目，k 点密度决定布里渊区采样的分辨率，这两个参数互相独立，但共同影响计算精度。先测截断能，再测 k 点，是大多数从业者固定的顺序：固定一个较大的 k 点网格，把截断能从 400 eV 加到 600 eV，观察总能量变化，取能量差小于 1 meV/atom 的最小值。然后固定这个截断能，增加 k 点密度，直到总能变化同样低于阈值。

这里有个取舍原则：截断能影响的是单个原子附近波函数描述，k 点密度影响的是周期性系统的电子分布。不同体系敏感度不一样。比如氧原子的小核赝势对截断能很敏感，而空旷的层状结构对 k 点的 c 方向密度不敏感。用工具箱的好处是，扫描参数时只需要改一个参数值，重新生成输入即可，不需要动其他文件。

带隙值对 k 点密度尤其敏感，直接使用默认网格可能高估带隙。对于带隙算不准的体系，把 k 点加密一档后带隙变化超过 0.05 eV，就要继续加密。如果想用杂化泛函，截断能也要相应提高，通常比 GGA 高出 30% 到 50%，否则交换积分部分描述不充分。

4.2 对称性开关与固定模式的取舍

对称性开关是个能省时间但也容易藏问题的选项。开启对称性后，程序只独立计算不可约布里渊区里的一部分 k 点，结果自动展开到整个布里渊区，计算量能降低一半甚至更多，所以默认是打开的。但开启对称性也有代价：如果输入的坐标精度不够，程序会把近似对称的结构误判成对称性较低的晶体，导致重复计算或者错误的对称性分类。

遇到这类情况我会在输入文件里把对称性开关显式设一下：结构来自实验数据且原子位置精度较高时开启对称性；结构来自猜测模型或原子位置经过随机扰动时，先关掉对称性跑一轮结构弛豫，再开启对称性算电子结构。这样一个流程能规避多数对称性误判问题。工具箱的输入函数也提供这个选项，默认开启，需要时手动关闭。

固定模式是另一个容易踩的选项。表面计算里，如果只固定中间几层原子让它驰豫，其余层保持理想构型，工具箱需要在结构体里额外标识哪些原子是活动的。这个标识用整数数组实现：值为 1 表示该原子参与弛豫，值为 0 表示固定。如果漏了这个标识，整个结构会被全部放松，表面模型在计算中逐渐变形，最后得到一堆没有意义的数据。我的习惯是先检查固定原子数量再提交任务。

4.3 批处理参数扫描与收敛判据

写脚本批量扫描参数是这个工具箱优势最大的场景。以下脚本对一组截断能逐一生成输入文件、提交计算、并读取总能量：

base = dfttoolbox.input.vasp(cryst, ... 'encut', 400, 'kpoints', [9 9 3], ... 'xc', 'PBE', 'sigma', 0.1, 'ismear', 1); encut_list = [400 500 600 700 800]; E_summary = zeros(length(encut_list), 2); for ie = 1:length(encut_list) folder = sprintf('scan_encut_%d', encut_list(ie)); inp = base.clone(); inp.set('encut', encut_list(ie)); inp.write(folder); run_job(folder); % 调用集群提交命令，按你的环境修改 E_summary(ie, :) = [encut_list(ie), ... dfttoolbox.result.read_final_energy(folder)]; end % 打印收敛结果，单位为 eV/atom E_ev_atom = E_summary(:, 2) / numel(cryst.species); fprintf('encut E/atom: \n'); disp([E_summary(:, 1), E_ev_atom]);

代码逻辑：先确定基准结构和其他参数，再循环不同截断能。clone 方法复制基础输入对象，set 方法覆盖单个参数，这样其他设置保持一致。run_job 需要你根据自己环境替换成实际的提交命令，本地跑就直接调用系统命令。

扫描完后用能量差判断收敛。比如 500 与 400 eV 差 12 meV/atom，600 与 500 差 1.5 meV/atom，700 与 600 差 0.3 meV/atom，取 600 eV 就足够，没必要用 800。同理对 k 点网格从 [7 7 3] 扫到 [11 11 5]，每个方向密度翻倍，直到总能差小于 1 meV/atom。注意扫描过程里不要同时调两个参数，否则无法定位瓶颈，白跑一堆作业。

这里的资源开销是可以预估的。一个 10 原子左右的体系，k 点从 [7 7 3] 加到 [11 11 5]，单次耗时大约翻一倍半。如果发现某个方向上加密后结果几乎不动，就说明已经收敛，可以停止扫描。参数扫描的最终产物是一个收敛参数表，之后所有同类型体系都用这组参数，能省下大量重复测试时间。

5. 避坑：DFT 结果分析最常踩的五个现场

5.1 费米能级漂移，能带整体错位

现象：能带图画出来导带底比预期低 0.3 eV，态密度里费米能级附近的峰位置和实验对不上。

原因：费米能级是从计算日志里读的，但读取的位置和能带文件里的能量零点不一致。某些程序能带文件默认以费米能级为零点，另一些则输出绝对能量，手动混用后所有能带整体偏移。还有一种情况是自洽计算设置了固定占据数，费米能级本身没有经过展宽处理，读出来的数值不是真正的费米能级。

解决：统一从同一个输出文件解析费米能级，并跟能带数据里零点的约定对齐。我一般会先打印几行原始能带数据，手动算一下价带顶位置，与工具箱读出的费米能级对比，确认零点一致后再批量处理。固定占据数的计算需要先跑一遍自洽拿到费米能级，再把它传给能带计算，直接拿上一次的费米能级复用不靠谱。

5.2 态密度曲线上的诡异尖峰

现象：态密度曲线在某个能量位置出现一根又窄又高的尖峰，能量间隔只有 0.01 eV，周围数据点却正常。

原因：输出的能量网格不均匀，读取时用了固定展宽的重采样，把稀疏区域的高值点放大。更常见的一个原因是 k 点网格太粗，导致某些 k 点上存在几乎平的能带，态密度在这些位置出现数值上的表观发散，物理上并不代表尖锐的电子态。

解决：先增加 k 点密度重新计算，确认尖峰是否随网格加密而变弱。如果尖峰消失，说明是采样不足的假象；如果尖峰依然存在，检查这个能量位置是否对应某条能带的范霍夫奇点。工具箱里可以把重采样展宽增大再看峰形，但调展宽只是掩盖问题，不是根治方案。

5.3 原子坐标重复导致对称性误判

现象：输入结构明明是对称的，但程序算出来的对称性降到最低，计算量暴增，结果还出现虚构的原子位移。

原因：结构文件里坐标精度只给到小数点后三位，晶格常数给到 0.01 Å，程序在做对称性判定时发现原子间距差超过容差，直接放弃对称运算。还有一种情况是分数坐标恰好越过周期性边界，同一原子在文件里出现了两份坐标，导致原子数翻倍。

解决：把原子坐标精度统一到 6 位小数以上，晶格常数保留 3 到 4 位，然后重新生成输入。对于周期性边界的问题，按最近镜像距离检查一遍坐标，确认所有原子都落在 0 到 1 的区间。工具箱读取结构时会做坐标归一化检查，但外部来的结构文件还是需要在写入前用可视化程序核对一遍原子位置。

5.4 总能量一直在降却没有收敛

现象：能量收敛曲线显示每个离子步的总能量都在降低，但跑了 40 步还达不到收敛标准，电子步和离子步数目也异常多。

原因：最常见的是初始位移设置不当，结构弛豫在极小值附近来回震荡，每个离子步只下降一点点。另一种原因是截断能过低，原子移动时能量变化非常平缓，力收敛缓慢。

解决：把离子步最大步数调小，先跑到 20 步左右，观察每一步的位移是否在减小。若位移降到接近阈值又反弹，说明初始构型离平衡太远，可以先做一次粗弛豫，再拿弛豫后的结构继续精修。截断能过低的情况，按第 4 章的流程先测截断能，不要盲目增加离子步数。血泪教训是曾经把步数设到很大，多跑了两天，最后发现问题只是截断能低了一百多电子伏特。

5.5 把全部波函数塞进 cell 数组导致内存爆掉

现象：读大体系的波函数或本征矢量时，MATLAB 直接报内存不足，或者反复触发交换，整台机器卡死。

原因：波函数矩阵大小是 k 点数 × 带数 × 平面波基组数，大体系下轻松超过几个 GB。工具箱读取时默认只保留每个 k 点的能量本征值，但如果读取函数开启了输出全部波函数系数，就会一次性把所有数据加载到内存。

解决：按需读取，只保留能带能量和态密度数据，不要勾选波函数输出选项。如果确实需要波函数做投影，用分块读取的方式，逐 k 点处理之后再释放变量。MATLAB 里用 clear 显式删除临时变量，或者把读波函数的操作包在独立函数里，让局部变量随函数退出自动释放。内存问题的第一道防线是搞清楚数据量级再决定要不要读。

6. 进阶：把工具箱封装成批量流水线，跑完一组结构自动出对比报告

6.1 批量提交与结果收集框架

当要比较多个结构的稳定性时，逐个手动跑显然不现实。我会把前面的输入生成、参数扫描和结果读取串成一条流水线：一个目录对应一个结构，循环完成输入构建、任务提交、结果收集与保存。核心代码只有十几行，但每天能省出大块时间。

materials = {struct('a', 3.85, 'b', 3.85, 'c', 9.42, ... 'species', {{'Ti', 'O', 'O'}}, ... 'coords', [0.3333, 0.6667, 0.25; ... 0.6667, 0.3333, 0.75; ... 0.5, 0.5, 0.5]), ... struct('a', 4.02, 'b', 4.02, 'c', 10.60, ... 'species', {{'Zr', 'O', 'O'}}, ... 'coords', [0.25, 0.75, 0.5; ... 0.75, 0.25, 0.75; ... 0.5, 0.5, 0.25])}; for i = 1:numel(materials) folder = sprintf('material_%02d', i); inp = dfttoolbox.input.vasp(materials{i}, ... 'encut', 520, 'kpoints', [9 9 3], 'xc', 'PBE'); inp.write(folder); run_job(folder); result{i} = dfttoolbox.result.collect(folder); save(sprintf('result_%02d.mat', i), 'result'); end

这里的 cell 数组 materials 存放多个结构体，循环里每个结构独立建目录，避免文件相互覆盖。save 每一步都做一次快照，即使中途任务失败，也能从最近一次保存点接着分析，相当于留了后悔药。

6.2 一键导出对比图与数据表

收集完所有结果后，我会把总能量、带隙和态密度峰位整理成一个表格，并用工具箱的汇总函数直接导出图和数据文件。这样拿着一个脚本就能给出格式统一的结果，不用在报告前临时整理一堆散落的曲线。

summary_table = dfttoolbox.report.summary(result, ... 'properties', {'energy', 'gap', 'dos_peak'}, ... 'outfile', 'summary.xlsx'); dfttoolbox.report.figure_energy_comparison(result, 'energy.png');

summary 函数会自动提取每个任务的收敛能量、带隙和态密度峰值，写成表格；figure_energy_comparison 把多个结构的能量做横向对比，适合直接放进报告。

6.3 一个让我长了记性的习惯

写完脚本先拿一个小体系做全流程验证，确认输出文件格式、单位、数据维度无误，再运行大批量任务。这个习惯让我避免了一次把 20 个结构提交到集群后发现输入少了 k 点设置的失误。从那以后我每次跑 DFT 任务前，都强制走一遍这段脚本化流程，确认每一步都能生成可复现的结果，希望帮到你。

本文还有配套的精品资源，点击获取