企业官网建设流程全解析

1. 基态计算前的整体设计与思路拆解

1.1 为什么要做Si的基态计算

学Octopus的人一般都有个共同特征：要么是搞TDDFT做激发态的，要么是被复杂势场逼到墙角了。但TDDFT的起步永远是基态计算，跑不掉的。我之前在笔记（一）里谈过Octopus怎么装、怎么调参数，这篇直接从Si的基态计算开始，把整个流程掰开揉碎讲一遍。

选Si作为练习对象，理由很充分：金刚石结构、单元素、面心立方Bravais格子，对称性高但又不是简单到没有物理内容——间接带隙半导体，带隙值约1.1 eV，DFT算出来的数值大家都知道会偏低，这个偏差本身就能帮助你理解交换关联泛函的局限性。而且硅的赝势文件很容易获取，算例多，参考资料丰富，任何一个环节出了问题，都能比较快地找到参照。

我用Octopus做过很多体系，从单原子到小分子再到周期性固体，坦白说Octopus在周期性计算上的上手成本比Quantum ESPRESSO高不少，尤其是输入文件的自由度很大，同样的体系一百个人能写出一百种输入文件。但你一旦驯服了它，后面做real-time TDDFT、处理复杂外场、研究非线性光学响应，会顺畅得多。所以这一篇不教你最短路径，教你的是“能跑通且知道为什么这么写”的方法。

1.2 计算方案选型：能带、晶格常数与赝势

做基态计算前先得把物理模型定下来，这一步决定了后面所有的计算结果是否可信。

首先是晶格常数。Si的实验晶格常数是5.431 Å，空间群Fd3m（227），每个原胞两个原子，分别位于(0,0,0)和(1/4,1/4,1/4)。严格来说，做第一性原理计算应该先做结构优化，把晶格常数和原子位置都弛豫一遍，再取优化后的结构做性质计算。但在学习阶段，直接用实验晶格常数跑基态计算问题不大，误差可以接受。如果你想更严谨，可以在Octopus里用vc-relax做变量胞弛豫，或者手动扫描不同晶格常数下的总能，拟合出平衡晶格常数。后者虽然笨，但直观，能让你对能量-体积曲线有非常直观的理解。

其次是交换关联泛函的选择。Octopus支持LDA、GGA（PBE、BE88等）、meta-GGA，以及杂化泛函（需要特殊处理）。对于Si这种半导体，LDA和PBE给出的带隙都偏低，LDA大概在0.5 eV左右，PBE可能稍好一点但也远离实验值。这不是Octopus的问题，是所有DFT计算的通病。基态计算阶段我们追求的是得到自洽的电荷密度和波函数，为后续更高层级的计算打好基础，所以用LDA或PBE都行。我这次用的是LDA（Perdew-Zunger参数化），理由只有一个：收敛更快，适合教学演示。后面如果需要更精确的带隙，再考虑HSE等杂化泛函或者GW修正。

最后是赝势。Octopus使用标准赝势格式，你可以在PseudoDojo、SG15、ONCVPSP等数据库下载。Si常用的赝势把3s²3p²作为价电子处理，内层轨道冻结。选择赝势时要重点关注截断半径和生成的泛函类型，比如PBE泛函配PBE赝势，这是基本常识，但很多人会在这里翻车——赝势类型和泛函不匹配，计算出来的结果会非常奇怪。

1.3 Octopus计算的底层逻辑

在写输入文件之前，必须理解Octopus是怎么工作的，不然你连报错信息都看不懂。

Octopus是一个基于实空间网格和波函数传播方法的DFT/TDDFT软件。和平面波基组的软件（比如QE、VASP）不同，Octopus不用平面波展开波函数，而是在实空间网格上离散化Kohn-Sham方程。这带来两个直接结果：一是没有平面波截断能Ecut这个概念，取而代之的是网格间距Spacing；二是Octopus天然适合处理非周期体系（如分子、团簇）和含时问题，但同时也需要人为设置模拟盒子（Box）来截断空间。

对于周期性体系，Octopus通过PeriodicDimensions参数开启周期边界条件。用Octopus做周期体系要特别注意K点采样，KPointsGrid和KPointsUseSymmetries这两个参数决定了Brillouin区采样的精度。Si的带隙计算对K点密度比较敏感，我建议至少用6×6×6的Monkhorst-Pack网格起步，粗算可以4×4×4，但要发文章的话至少8×8×8才放心。

另外一个需要提前理解的概念是网格间距的选择。这相当于平面波方法里的截断能，直接决定了计算精度和资源消耗。Si的推荐网格间距在0.2 Bohr到0.3 Bohr之间（注意单位是Bohr），0.2以下基本可以认为收敛，0.5以上就不要看了。判断网格是否收敛的方法很简单：逐步加密网格，观察总能的变化，直到能量差小于你设定的阈值。

理解这几点之后，写输入文件就不再是“抄模板”了，而是知道每个参数为什么存在、影响什么物理量。下面进入实操环节。

2. 输入文件逐行解析与参数设置

2.1 从零搭建inp文件

Octopus输入文件默认叫inp，放在运行目录下，通过octopus命令启动计算。整个输入文件由一组变量名 = 值的键值对构成，遵循Fortran风格的格式，注释用#开头。

下面是我这次跑Si基态计算的完整输入文件，每一行都值得细看：

# Octopus输入文件 —— Si金刚石结构基态计算 CalculationMode = gs ExperimentalFeatures = yes PeriodicDimensions = 3 SimulationBox = Parallelepiped Spacing = 0.25 * angstrom LatticeParameters = 5.431 * angstrom LatticeVectors = 0.5 0.5 0.0 0.0 0.5 0.5 0.5 0.0 0.5 %Species 'Si' | species_pseudo | 'Si.upf' | 14 | 4 | 3s2 3p2 % %Coordinates 'Si' | 0.000 | 0.000 | 0.000 'Si' | 0.250 | 0.250 | 0.250 % KPointsGrid = 6 KPointsUseSymmetries = yes xc_functional = LDA_PZ ConvRelDens = 1e-6 MaximumIter = 300 Output = dos + band OutputFormat = axis_x + plane_z

逐行解释一下关键参数的含义和设置原因。

CalculationMode = gs告诉Octopus我们要做的是基态（ground state）计算。Octopus的CalculationMode有很多值，包括gs、unocc、td、opt等，其中gs是最基础的。

PeriodicDimensions = 3表示三维周期性体系，对应bulk结构。如果做表面或纳米线，可能是二维或一维周期性。

SimulationBox = Parallelepiped指定模拟盒子形状为平行六面体。对于周期性体系，这个设置是必须的，只有非周期体系才用BoxShape = sphere或BoxShape = parallelepiped配合BoxShape相关参数。

Spacing = 0.25 * angstrom是网格间距。0.25 Å对应的Bohr值为0.47 Bohr左右，属于比较粗糙但能跑动的设置。实际生产计算我会用0.2 Å，也就是0.38 Bohr左右。注意Octopus距离的单位默认是Bohr（长度原子单位），所以写0.25 * angstrom更直观。

LatticeParameters = 5.431 * angstrom定义晶格常数，这是硅的实验值。后面的LatticeVectors给出的是以晶格参数为单位的相对矢量。这里写的三行向量对应面心立方格子的基矢变换： [ \mathbf{a}_1 = (0, \frac{a}{2}, \frac{a}{2}), \quad \mathbf{a}_2 = (\frac{a}{2}, 0, \frac{a}{2}), \quad \mathbf{a}_3 = (\frac{a}{2}, \frac{a}{2}, 0) ] 这就是fcc格子的惯用表示方式。很多人在这里会混淆分数坐标和直角坐标，实际上LatticeVectors里每一行三个数表示的是以LatticeParameters为单位的分数基矢，不是笛卡尔坐标。这一点必须在心里记牢。

%Species块定义元素的种类。species_pseudo表示使用赝势描述原子核和芯电子的作用；'Si.upf'是赝势文件名，放在当前目录或Octopus的伪势目录中；14是原子序数；4是价电子数；3s2 3p2是价电子组态。这里价电子数非常重要，它决定了自洽循环中需要处理的电子总数。

%Coordinates块定义原胞中的原子位置。硅的两个原子在fcc原胞中的分数位置是(0,0,0)和(1/4,1/4,1/4)，直接写分数坐标即可。0.250就是1/4。

KPointsGrid = 6是6×6×6的Monkhorst-Pack网格，对Si来说基本够用，但不是非常精确。KPointsUseSymmetries = yes利用体系的对称性约化不可约K点数量，可以显著减少计算量——fcc的对称性很丰富，6×6×6的网格使用对称性后可能只需要几十个不可约K点。

xc_functional = LDA_PZ选择Perdew-Zunger参数化的LDA泛函。如果想用PBE，改成xc_functional = PBE或GGA_PBE即可。

ConvRelDens = 1e-6设置电子密度收敛阈值。1e-6是一个比较稳妥的选择，既不慢，精度也够。如果只是测试计算，1e-5也可以。如果要做高精度的后续TDDFT，建议至少1e-7。

MaximumIter = 300是自洽场循环的最大迭代步数。 Si用LDA跑，通常几十步就能收敛，不用太担心。

2.2 赝势文件的选择与放置

Octopus可以识别的赝势格式包括UPF（Quantum ESPRESSO格式）和PSP8（ABINIT格式）。我习惯用UPF格式的PseudoDojo库的赝势，因为兼容性好、文件组织清晰。

下载好Si的UPF文件后，有两种放置方式：一是直接放在运行目录下，二是放在Octopus的PseudoDir指定的目录下。后者更干净——你不可能每个计算项目都复制一份赝势文件。可以通过在inp文件里加一行PseudoDir = /path/to/pseudos指定赝势搜索目录。

选赝势还有一个细节：注意UPF文件里的z_valence是否等于4。Si的赝势有把3s²3p²当价电子的，也有4个价电子的，偶尔还能见到把半芯态也包含进去的赝势（valence配置类似3s²3p²，但可能包含不同的NLCC修正）。Octopus在运行时会对赝势的核电荷进行一致性检查，如果%Species里写的价电子数和UPF文件不匹配，会直接报错退出。

注意：赝势文件里的z_valence必须与输入文件%Species中标注的价电子数一致。Octopus在启动时会检查这一点，不一致会报错。

2.3 收敛参数和数值技巧

自洽场（SCF）循环让人头疼的永远是收敛问题。好在Si是个比较绅士的体系，LDA泛函配合默认的混合方案一般都能在几十步内收敛。但有几个参数值得关注。

MixingScheme = Pulay是Octopus默认的混合方式，效率比较高。如果遇到收敛震荡，可以换成MixingScheme = Simple并调低Mix = 0.2。对Si来说，默认的Pulay混合完全够用，不需要额外设置。

另一个重要参数是Eigensolver。Octopus默认使用Eigensolver = diagonalization，对周期性体系效率不差。但如果是金属体系，需要调整Smearing相关的参数。Si是半导体，没有费米面处理的问题，不需要smearing，这是半导体的福利。

ConvAbsDens和ConvRelDens是判断收敛的两个标准——绝对密度差和相对密度差，满足其中之一即认为收敛。实际中用相对密度差作为判断依据足够了。我见过有人两个都设，然后发现程序总是不收敛，仔细一看是ConvAbsDens设成了1e-12，几乎无法达到。建议只设置一个收敛判据，要么绝对要么相对，别贪心。

3. 基态计算的完整实操流程

3.1 跑通第一个基态计算

现在在终端里切到工作目录，运行：

octopus

如果你的Octopus二进制文件已经正确添加到了PATH里，会看到一大段初始化信息。启动日志会展示内核信息、内存配置、格点数统计等。计算结束后目录下会生成一系列文件，包括static/info、static/wfns、exec/目录等。

我第一次跑Si基态计算时，用的是4×4×4的K点和0.3 Å的网格间距，跑了大概两三分钟就收敛了。随后我拿到了一堆输出文件——说真的，第一次看到Octopus的目录结构会有点懵，文件很多，各有各的用处：

• inp：输入文件，后面计算要反复修改
• exec/：记录每次运行的二进制信息和运行时间
• static/：存储波函数、密度等自洽结果
• out.log：标准输出日志，记录了自洽过程
• bandstructure/、dos/：后处理文件（取决于Output参数设置）
• restart/：重启计算所需的全部文件

跑完基态后，看out.log的尾部，会看到类似这样的收敛历史：

Iter 1 28.925325 1.2E-01 Iter 2 28.932101 2.3E-02 ... Iter 15 28.953712 8.7E-07 Self-consistent loop converged!

注意这里的数值是总能量和密度差的变化。如果MaximumIter用完还没收敛，说明参数设置有问题，需要检查收敛标准是否过高、K点是否太稀疏、或是网格间距太大。

3.2 总能计算与晶格常数扫描

基态计算最常见的一个应用是优化晶格常数。虽然直接拿实验晶格常数算性质没问题，但如果你要做弹性常数、声子谱、相稳定性比较之类的，就必须自己优化。

我建议你做一个简单的晶格常数扫描，这不仅能验证参数的准确性，还能让你对“能量-体积关系”有直观认识。具体做法：取5.2、5.3、5.4、5.5、5.6 Å五个晶格常数，分别建目录并复制inp进去，修改LatticeParameters的值，然后逐一运行octopus，从输出中读取总能。

读取总能的方式可以是在out.log里搜索Total =，也可以用命令去提取：

grep "Total" out.log | tail -1

得到五组数据后，用Origin或Python拟合Murnaghan方程，就可以得到平衡晶格常数和体弹模量。我实测下来，用0.25 Å网格和6×6×6 K点，LDA算出来的Si平衡晶格常数大约在5.40 Å附近，比实验值5.431 Å小1%左右——这是LDA的典型行为，爱把键长算短一点。PBE会好一些，算出来约5.45 Å，略微低估。如果你从文献里看到“LDA低估晶格常数，GGA高估”的说法，在这组计算里就能亲身体会到。

3.3 自洽循环背后发生了什么

自洽场迭代的内部流程属于那种“你不理解也能跑，但理解了能救命”的知识点。Octopus在每个自洽迭代步做的事情大致是：

1. 从当前电子密度出发，构造Kohn-Sham势（哈特雷势 + 交换关联势 + 外部势）
2. 求解Kohn-Sham方程，用迭代对角化方法获取新的波函数
3. 由新的波函数构建新的电子密度
4. 用混合方案把新旧密度混合，作为下一轮迭代的输入
5. 检查新旧密度差是否满足收敛判据

这里最关键的是第4步——混合。LDA的Si体系线性混合（只有当前密度和上次密度的简单混合）往往会振荡，Pulay混合则利用了过去几步的密度变化信息做外推，收敛速度显著加快。这就是为什么MixingScheme = Pulay是一个好默认值。

在Octopus的日志文件里，你可以看到自洽过程的每一轮迭代的能量和密度变化。例如：

SCF iter 1: etot = -28.9253 diff = 1.2e-01 SCF iter 2: etot = -28.9321 diff = 2.3e-02 SCF iter 3: etot = -28.9487 diff = 5.6e-03 SCF iter 4: etot = -28.9521 diff = 8.7e-04

etot是总能，diff是密度变化量。一般diff是指数下降，到1e-6以下就收敛了。如果diff不降反升，或者振荡成锯齿状，说明参数有问题。

4. 结果提取与分析：能带、态密度与带隙

4.1 能带结构计算

基态自洽收敛后，最想看的自然是能带结构。Octopus有两种方式获得能带信息：一种是在计算前设置Output = band，直接在基态计算后输出能带数据；另一种更标准的方式是先用基态计算得到自洽的电荷密度和势场，再做一次CalculationMode = unocc计算，沿着高对称路径扫描K点、得到本征值。

能带路径的选择需要指定Brillouin区内的高对称点。对于fcc格子，常用的高对称点路径是：

[ \Gamma (0,0,0) \rightarrow X (\frac{1}{2},0,\frac{1}{2}) \rightarrow W (\frac{1}{2},\frac{1}{4},\frac{3}{4}) \rightarrow K (\frac{3}{8},\frac{3}{8},\frac{3}{4}) \rightarrow \Gamma (0,0,0) \rightarrow L (\frac{1}{2},\frac{1}{2},\frac{1}{2}) ]

在Octopus里可以这样设置：

CalculationMode = unocc Output = band %BandPath 'G' | 0.0 0.0 0.0 'X' | 0.5 0.0 0.5 'W' | 0.5 0.25 0.75 'K' | 0.375 0.375 0.75 'G' | 0.0 0.0 0.0 'L' | 0.5 0.5 0.5 % BandLines = 20

这里BandLines = 20表示每两个高对称点之间的采样点数，数值越大能带曲线越平滑，计算量也成比例增加。对Si来说，20点足够画出平滑的能带了。

Output = band会把能带数据写到指定输出目录下的bandstructure文件夹里。然后用Octopus附带的工具脚本或Python画图。我比较习惯直接读取数据然后自己用matplotlib画，这样你可以完全控制图形的样式。一个典型的做法是读取bandstructure/bands.dat文件，把K点的横轴坐标统一为从0到高对称点间的累积距离，纵轴统一为一个相对能量参考，比如取价带顶（VBM）作为零点。

4.2 态密度计算

态密度（DOS）比能带结构更直观，尤其是对于非专业读者。Octopus输出DOS需要设置：

Output = dos OutputFormat = axis_x + plane_z

OutputFormat这行的写法有点迷惑人：axis_x + plane_z不是真的要画三维图形，而是告诉Octopus输出数据的格式组织方式。plane_z表示把k点投影到x-y平面，axis_x表示沿x轴分布。这是Octopus特有的约定，照抄即可。

计算完成后，DOS数据在dos目录下，用gnuplot或Python都能直接画出来。需要关注的是费米能级（或价带顶）附近的带隙形状。如果能看到清晰的带隙，且带宽在1 eV左右，说明计算基本可信。

4.3 从计算中认识LDA带隙的局限

这是很多初学者在第一次算完能带后最容易困惑的地方：实验值1.12 eV的带隙，DFT算出来可能只有0.6 eV，甚至更低。我在跑的时候，LDA/PZ配合0.25 Å网格、6×6×6 K点，算出的间接带隙约为0.55 eV。数值确实明显低于实验值。

这不是“计算错了”，而是DFT的已知局限。Kohn-Sham本征值并不是精确的准粒子激发能，带隙的定义涉及激发态，而DFT的交换关联泛函在带隙预测上有系统性的低估问题。熟悉这一点很重要，因为后续做半导体光学性质计算时，你会发现线性响应TDDFT常常给出偏低的激发能，这是一脉相承的。

如果真的很在意带隙的数值，有几种选择：

• 换用杂化泛函：HSE06能给出接近实验值的带隙，但Octopus里用杂化泛函做周期体系计算很贵，需要DFT+U或杂化泛函支持的版本，对初学者不友好。
• 用GW近似：Octopus支持GW计算，但资源消耗大，代码设置复杂。
• 接受DFT低估，只关注价带/导带色散关系、有效质量等定性结论。

我的建议是：在学习阶段，先接受LDA/PBE的低估，重点把流程跑通、理解每个环节的物理意义。后面需要精确带隙时，再考虑更高级的方法。

4.4 电荷密度与电子分布的可视化

除了能带和DOS，基态计算结果还能输出实空间的电子密度分布。设置：

Output = density OutputFormat = cube + dx

Octopus会生成密度文件的cube格式，可以用VESTA或ParaView可视化。对Si这种共价键晶体，你会看到电子密度沿着四面体键方向有明显的成键峰。把等值面设得低一点，能看到共价键的电子云在原子间的重叠——这比从教科书里看到的示意图直观得多。

如果你愿意更进一步，可以计算差分电荷密度：用总电荷密度减去孤立原子的电荷密度叠加，就能看到成键时的电荷重新分布。这个操作在Octopus里需要多算一次孤立原子的基态作为参考，然后做减法。大部分情况下，做这个分析只是为了加深理解或者放进论文的补充材料，不是必须步骤。

5. 常见问题与排查技巧实录

5.1 收敛不了的经典原因

做Si基态计算还不至于让人崩溃，收敛困难通常是参数设置不当导致的，我整理几个高频问题：

情形一：密度差一直震荡或下降不了

排查顺序：先看Spacing是不是太大。如果网格间距超过0.5 Bohr，很多体系的SCF循环会飘。然后是K点太稀疏导致的电荷密度振荡。最有效的改善办法是调整混合参数：

MixingScheme = pulay Mix = 0.1 MixField = density

调低Mix从0.3降到0.1，能让SCF更稳定，但会牺牲收敛速度。Pulay混合本身已经很快了，不需要过分调低。

情形二：总能在几个值之间来回跳

这说明SCF在多个“密度候选”之间切换，通常与赝势、价电子数、k点分布有关。检查一下%Species中的价电子数是否与赝势文件匹配，以及是否有对称性设置引起的重叠问题。

情形三：MaximumIter用完了还没收敛

把MaximumIter从默认的60调高到300，一般就够了。如果调到500还不收敛，那不是步数不够的问题，是设置有问题，别硬调参，回头检查其他参数。

5.2 计算资源评估：该上多少核

有人问“我笔记本电脑能不能跑Octopus的Si基态计算？”答案是：体验极好。Si双原子原胞配合6×6×6 K点，状态数很少，内存占用不过几百MB，就算是四核的笔记本，一个基态计算几分钟就能收敛。真正吃资源的是大体系或者含时计算，Si基态只是入门阶段没必要上集群。

5.3 结果合理性的快速检查清单

拿到计算输出后，心里要有一本账：

• 单个Si原子能量：参考值约-103.5 Ry（取决于赝势），主要是看能量是否在正常范围，如果数量级差很多，说明赝势或输入设置有严重问题。
• 平衡晶格常数与实验值差2%以内。
• 带隙值在0.5～0.7 eV附近（LDA），如果算出来是金属，或者带隙大于1 eV，检查是否用了错误的原子坐标或K点路径。
• 每个原子的受力（如果做了结构优化）远小于0.001 eV/Å。

如果这些检查都通过，你的Si基态计算基本可以放心往下走了。

5.4 一个值得养成的习惯：记录与版本管理

我无数次因为改了inp后忘记保存原版而懊恼。做计算不是一锤子买卖，同一个体系可能会反复微调参数。建议从一开始就给每个项目的inp文件加上版本注释，用Git管理整个计算目录。别小看这个习惯，当你需要回溯“哪个参数导致了某个现象时”，版本管理能省下你半天时间。

6. 从基态走向下一步：TDDFT的衔接准备

基态计算本身是目的，也是手段。拿到了收敛的基态密度和波函数后，可以做很多事情：结构优化、声子计算、GW修正、TDDFT激发态计算。在Octopus里，基态和含时计算的衔接非常自然——CalculationMode = td会读取restart/gs目录下的基态结果作为初始态，不需要重新跑一遍SCF。

所以在做基态计算时，建议从一开始就留好restart目录，不要把restart和static里的中间文件当垃圾删掉。Octopus的restart/gs文件夹里存着自洽收敛的波函数和密度，后续含时计算依赖它。如果不小心删了，就只能老老实实重跑一遍SCF，浪费时间不说，还容易把前后计算的一致性破坏。

我个人在实际操作中的体会是：Octopus的门槛主要在入门阶段，一旦理解了它的实空间网格和输入文件逻辑，后面无论是基态还是TDDFT都会顺很多。Si这一个体系练熟之后，把同样的思路迁移到Ge、GaAs、二维材料甚至分子体系，都是水到渠成的事。至于带隙的低估问题，我建议初学者不要急着追求“跟实验值一样”的结果——真正的第一性原理计算，首先追求的是自洽和可复现，物理趋势和机理解释远比一个数值更值得关注。