“NI”这个符号一出现，搞计算材料的基本都明白是什么意思。镍，周期表第28号元素，3d过渡金属里的老常客，催化、磁性、电池、合金哪哪儿都有它。但真跑到DFT平面波代码里算镍的时候，事情就没那么舒坦了——3d电子局域性强，和平面波基底天生不对付，你要是随便拿一套赝势糊弄过去，后面无论是算晶格常数、磁矩还是表面吸附能，都会给你埋一堆雷。
ATOMPAW生成PAW_PBE赝势这件事，本质上是把“原子尺度的全电子量子力学问题”和“晶体尺度的平面波计算”做一个干净利落的接口。PAW方法，Projector Augmented Wave，投影缀加波，是现在平面波DFT代码里最主流的方案之一。而ATOMPAW就是那套能从原子全电子计算出发，手搓出PAW赝势数据的开源工具链。它不依赖哪个商业代码的闭源赝势库，生成过程透明，参数自己可控，最重要的是，你可以按需生成PAW_PBE赝势——也就是配PBE交换关联泛函的那一套。
这篇文章我就用Ni当靶子，完整盘一遍从电子构型分析、原子数据准备、截断半径选择，到ATOMPAW实际运行、赝势质量检查和排错的全过程。适合正在折腾VASP、Quantum ESPRESSO、ABINIT这类平面波代码，又对官方赝势库的“黑盒参数”不放心的人。哪怕你之前没碰过ATOMPAW，只要知道一点DFT基础，按着后面的思路走一遍，也能搞出靠谱的Ni赝势来。
1. 整体设计与核心思路：为什么拿Ni开刀最合适
1.1 镍的电子结构难点在哪
镍的基态电子排布写出来是 [Ar] 3d8 4s2，但固体镍里实际上3d和4s轨道之间的电荷转移很敏感，经常出现3d9 4s1的“有效”构型。也就是说，你在原子阶段选什么参考构型、初态构型，都会直接影响赝势在固体环境中的行为。
更深一层的问题是，3d轨道在径向尺度上比4s要紧凑得多，空间上又和3s、3p半芯态部分重叠。平面波展开要求所有价电子波函数在截断能范围内足够平滑，可3d轨道那几根尖峰，必须在PAW框架里靠原子轨道投影来恢复。这就是为什么对Ni来说，传统的模守恒赝势往往需要很高的截断能（容易超过600-800 eV），而PAW方法可以在这个问题上松一大口气。
另外，Ni是有磁性的。铁磁态下Ni的磁矩大概是每原子0.6 μB左右，这个磁矩对电子结构极其敏感。你用不同的赝势生成方案，哪怕只是截断半径差0.1 bohr、或者芯态/价态划分差一个3p轨道，最后的磁矩都可能差出0.05-0.1 μB。对很多表面催化计算来说，这个差异直接决定吸附能能不能比得准实验。
1.2 PAW相比模守恒与超软赝势的机制优势
先快速梳理一下PAW的底层逻辑，因为ATOMPAW的所有参数设置都是围绕这套逻辑展开的。
模守恒赝势的思路是直接把原子核和芯态一起抹掉，用一个平滑的赝势替代，让价电子波函数在外面和真实波函数一模一样。但在处理3d这种强局域轨道时，为了达到“一模一样”，截断能往往被推得很高，计算成本爆炸。
超软赝势（USPP）走的是另一条路：允许赝波函数在芯区不守恒，用多一套补偿电荷来修正，代价是引入额外的投影算符和电子密度修正项。PAW则可以理解成“把全电子波函数找回来的超软赝势”——它保留了一个原子轨道基组展开的变换关系，真正把全电子波函数显式重建出来。
对Ni这种体系，PAW的好处特别明显。一个原因是3d态的化学环境变化很大，从金属到氧化物到表面吸附体系，d带的位置、占据都在变，PAW通过原子投影通道的灵活扩展，能更好保持“可转移性”。另一个原因是PAW计算出的磁矩和全电子结果更贴合，超软赝势在磁矩上经常出现莫名其妙的偏离，PAW就好得多。
用一句话给不熟悉的人解释的话：模守恒赝势是“换了一个假的势，让它产生和真的一样的波函数”，超软是“换了个假的波函数，再用补偿电荷找补”，PAW是“波函数暂时换成好算的，算出结果后把真实的芯区细节装回去”。
1.3 为什么固定选PBE泛函
项目的名字里带着PAW_PBE，那泛函就是Perdew-Burke-Ernzerhof的GGA泛函。这几乎是现在凝聚态DFT计算的默认选项，没有太多争议空间。
PBE对金属键、表面吸附、大部分氧化物结构都能给出一个和实验比较接近的结果。对于Ni，无论是晶格常数、弹性模量还是磁矩，PBE通常给得都比LDA更稳。LDA普遍高估结合能，把晶格常数压得偏小，Ni这种磁性金属更是重灾区。PBE虽然也略微高估晶格常数，但误差范围在1%以内的情况非常多见，文献可比性也强。
还有一个非常实际的原因：VASP自带的PAW赝势里所谓“PAW_PBE”其实就是用类似方法配合PBE生成的，Quantum ESPRESSO的pslibrary里也是PBE占大头。所以你手搓一套PAW_PBE赝势出来，可以直接用别人的文献参数对比验证，不至于因为泛函不一致而鸡同鸭讲。
ATOMPAW本身也内置了PBE的处理路径，你只要在模版里选对泛函，生成的径向数据和交换关联势都是配套好的，不用自己额外去查交换关联的数值表。
2. 生成Ni赝势前必须想清楚的四个物理参数
2.1 参考态与初态怎么定
PAW赝势生成的第一步，是给原子作全电子计算。这里需要两套电子构型：
- 参考态配置（reference configuration）：被用来构建PAW基函数和投影算符的那个原子态，相当于“锚点”。
- 初态配置（initial configuration）：生成电子密度初始猜测时用的，通常选一个容易收敛的离子态。
对Ni，常见的选择是参考态3d8 4s2，初态3d9 4s1。为什么这样挑？因为参考态要尽量接近晶体里的实际电子分布。固体Ni是d带接近填满、s带几乎退化的状态，3d8 4s2在各个体系间是比较好折中的一个选择。初态选3d9 4s1则是因为全电子计算里这个配置比较容易自洽收敛，特别在使用对数网格快速估算能级时，3d9 4s1比3d8 4s2的更少出现轨道翻转。
有些人会把参考态写成3d10，理由是Ni在化合物里经常出现d10闭壳层状态。但这样做会让4s态在原子层面的行为偏软，到金属Ni里容易出现微弱的不对劲。我自己的经验是，除非你专门针对Ni(II)配合物做高精度计算，否则3d8 4s2作为参考态是最省心的。
2.2 芯态和半芯态的取舍
在PAW生成里，“哪些电子算芯态、哪些算价态”这个问题，比很多人想象的更重要。
Ni的原子序数是28，1s到2p一共10个电子，这些离3d太远，肯定算芯态，没悬念。真正要掰扯的是3s、3p这8个电子。3s2 3p6的能量深度大约在-80到-100 eV量级，不算深，空间上和3d有相当大的重叠。
如果把这8个电子直接扔进芯态，计算会明显变软，平面波截断能要求也低，但后果是：镍原子在压力下的行为、晶格动力学、以及一些磁性性质会出问题。因为3p轨道和3d轨道的径向重叠在压缩环境下会被“压”出不可忽略的电子云重排效应，你要是把它们全当芯，等于假设这套电子云永远僵硬不动。很多标题里带“PAW_PBE Ni”的赝势在VASP库里有几种版本，其中常区分“3p作为价态”和“3p作为半芯态”的版本，根本原因就在这。
ATOMPAW里处理半芯态不算复杂，就是把3s、3p也一并列入价态集合，但要留意它们的径向波函数在芯区还有幅度，这会加大径向网格密度要求，也让最终的赝势文件大一点。配3p半芯态的Ni赝势在截断能相同的情况下会比纯3d4s价态版本稍微慢一些，测量下来大概多5%-10%的计算耗时，但对磁性体系的准确度收益完全值得这个代价。
2.3 截断半径（rc）怎么选
截断半径大概是整个ATOMPAW输入参数里最微妙的一个。它决定了全电子波函数在多大半径内被投影波函数替代，也决定了赝势的“软硬”程度。
- rc取得小（比如1.4 bohr）：全电子波函数被保留的范围大，精度高，更“硬”，但平面波截断能要求很高。
- rc取得大（比如2.2 bohr）：赝势很软，平面波收敛快，但会牺牲一部分对局域化学环境的敏感度，表面吸附、磁性分裂这类计算容易出偏差。
Ni比较常用的rc范围在1.6-2.0 bohr之间。比如VASP标准PAW_PBE里的Ni，对应的等效截断能在250-300 eV附近，这差不多就是rc取到1.7-1.8 bohr的感觉。如果你做高精度表面催化，我会建议rc往1.6再靠一靠，配合包含3p半芯态，会把吸附能计算里的基组重叠误差压得更低。
另外径向网格也不是越多越好。ATOMPAW默认使用对数网格，起点一般在10^-5 bohr量级，终点要能覆盖到rc以外足够远。对Ni，网格点取300左右通常够了；如果你发现生成的径向波函数在近核区出现非物理振荡，优先怀疑的是网格点数不足，而不是rc错了。
2.4 投影通道与角动量通道的设置
PAW投影算符的数量直接决定了精确度上限。每个角动量通道（l）至少需要一个投影函数，但对于Ni的d通道，强烈建议至少配两个通道。
原因在3d散射势的复杂性。d轨道的形状、能量位置在金属和绝缘体中变化很大，一个投影函数只能覆盖一定的能量范围，两个或三个投影函数才能把散射行为更好地插值到不同化学环境。具体到ATOMPAW的参数表格里，通常会用“每通道投影数目”和“局域轨道数目”这类选项，Ni的设置可以参考：s通道1个或2个投影，p通道1-2个，d通道2-3个。
有人会贪多，把所有通道都加满投影函数。确实更精确，但PAW投影矩阵会变大，计算成本上升，而且有些代码对过多投影通道支持得不太好。一个实践经验是：先跑一遍原子自洽，看看每个投影通道的占据数是否都在合理范围内（0-1之间）；如果某个通道的占据数突然接近2或者出现怪异的负值，那就说明投影通道设置过度了。
3. 实战：ATOMPAW生成Ni_PBE赝势的完整流程
3.1 环境准备与编译要点
ATOMPAW的编译不算难，主要依赖Fortran编译器和LAPACK/BLAS库。常见组合是gfortran + OpenBLAS，Intel机器上也可以用ifort + MKL，后者在大数组运算时略快几秒，但日常用真感觉不出来。
典型的配置流程就是在源码目录里运行configure，指定FC和LAPACK路径，make。之后会生成可执行文件atompaw，以及后处理工具genpaw（用于把中间数据转换为特定代码格式）。
如果你只需要VASP格式的POTCAR，编译时最好带上VASP相关的支持开关；目标代码是Quantum ESPRESSO就选qe的输出格式；ABINIT自然就是abinit格式。这些输出格式在配置阶段明确指定，后面生成转换时会顺很多。实际建议是：全都编上，反正源码里就那几条开关，多花不了两分钟。
3.2 输入文件长什么样
ATOMPAW的输入文件在不同版本之间略有差异，但核心内容块是一致的。下面我放一个针对Ni的参考模板，里面的字段顺序建议以你本机版本的官方文档为准——我写这个是让你看懂每块在干什么，而不是让你无脑照抄。
# ATOMPAW 输入示例：Ni，参考配置 3d8 4s2 # 字段说明：原子名 ｜ Z ｜ 处理状态数 ｜ 截断半径 Ni 28 4 1.8 # 每个态一行：主量子数 n，角动量 l，占据数 f，能量初始猜测 3 2 8 0.0 4 1 2 0.0 3 1 6 0.0 3 0 2 0.0 # 投影与局域通道设定 # 此处略——视版本不同，以交互输入或后续标志段完成第一行Ni是标识；第二行28是原子序数；第三行的4代表你要显式处理四个价态（含半芯态），千万别把那10个芯态也算进去。第四行1.8是截断半径，单位bohr。
接下来的四行分别是3d8、4s2、3p6、3s2。很多版本里每个状态还需要初始能量，填0.0通常表示让程序自动在参考态附近找初始能级。你要是遇到收敛困难，再手动给一个大约的轨道能量。
至于投影通道数量的给定方式，每个版本处理得不一样，有的放在交互提问里，有的放在文件后段。我习惯先用一个保守默认：s通道1个投影、p通道1个、d通道2个。后续如果做d带中心类的高精度计算，再把d通道加到3个重新生成一遍，对比一次便知差别。
3.3 运行ATOMPAW并理解输出
运行方式通常就一条命令：
atompaw < atompaw.in或者按你版本支持的参数方式：
atompaw -i atompaw.in -o atompaw.out程序启动后，会先做全电子原子自洽，然后在此基础上构造PAW数据集。整个过程从几秒到几十秒不等，如果你在HPC集群的登录节点上跑也基本无压力。
运行结束后，你会看到一批输出文件。核心的是原子波函数数据文件，里面包含：径向网格、各态的对数导数、投影函数矩阵、补偿电荷参数、交换关联势数据等。这些文件有些是给后续转换用的中间文件，有些是最终可以放进代码的赝势文件。
别急着打开二进制文件看——没意义。你应该看的是atompaw输出日志里那几行摘要：总能、各态占据数、投影函数的范数以及径向函数的平滑度指标。只要这些指标没报警，文件就是可用的。
3.4 转换成VASP/QE/ABINIT格式
同一套原子数据和投影参数，要配合不同DFT代码，还需要做格式转换。
- 目标代码是VASP：用genpaw或对应脚本生成POTCAR。VASP要求POTCAR里有完整的原子数据头、径向网格、投影函数，以及PBE交换关联势信息。生成后你可以用
grep "VRHFIN"看一眼头信息，确认元素、版本号这些元数据没乱。 - 目标代码是Quantum ESPRESSO：转为UPF格式。QE的UPF文件本质是XML，里面包含同样的数据集，你可直接在
pseudopotentials目录里引用。 - 目标代码是ABINIT：通常生成pspcode的PAW格式文件，再配合abinit的输入文件使用。
转换工具的具体命令每个版本有区别，但不外乎输入中间文件名+输出文件名。转换完成后，强烈建议做一个最小计算测试（一个Ni原子或一个Ni晶胞都行），把截断能设到400 eV左右，看看总能量能否在几个电子伏的合理范围内收敛。
3.5 输入参数走向验证：一个可落地的操作顺序
实际工程里，我建议你按这个顺序操作，不要跳步：
- 先用参考态3d8 4s2跑一遍全电子原子计算，确认总能和文献参考一致（Ni原子的PBE总能大约在-2700到-2750 eV量级，不同轨道数设置略有差异）。
- 用默认rc=1.8、默认投影通道生成第一版赝势。
- 用这版赝势算一个Ni的fcc晶胞，对比实验晶格常数（约3.52 Å）和磁矩（约0.6 μB/原子）。如果偏差大，先怀疑截断能不够，再怀疑投影通道或rc。
- 记录不同截断能下的总能量变化，画一个收敛曲线，确认你用于生产计算的截断能落在平台区。
- 如果对计算结果精度仍不满意，回到ATOMAPAW，调整rc或投影通道数量重新生成。
这套流程的优点是把每一次参数改动的效果都暴露在可观察的物理量上，而不是在“赝势文件本身好不好看”这个层面打转。
4. 赝势质量怎么检验才算过关
4.1 原子级检验：对数导数与全电子对比
ATOMPAW生成的日志里会包含各角动量通道的对数导数曲线。这是赝势质量的第一道关卡：在rc附近和几个典型能量点，赝势的对数导数应该和全电子波函数的结果重合得非常好。如果曲线在某个能量点出现明显分叉，说明这个通道的投影函数不足以描述散射，需要增加投影或调整rc。
另一个直观检验是把PAW计算出来的原子总能和真正的全电子计算对比。ATOMPAW在生成过程中本身也是在作全电子计算的，所以它内部就有一组基准数据。后处理工具里一般也带原子测试功能，能直接输出差值。一般来说，1 meV/原子以内的差异是合格的。
4.2 晶体级检验：晶格常数、磁矩、体弹模量
原子级检验通过后，真正有说服力的是晶体结果。
拿fcc Ni来说，实验晶格常数约3.52 Å，铁磁状态下每个原子的磁矩约0.6 μB。用PBE泛函配合靠谱的PAW赝势，计算结果通常在3.50-3.54 Å之间，磁矩在0.6-0.65 μB之间。假如算出来晶格常数偏到3.48以下或3.56以上，说明赝势的径向行为有偏，优先检查是否漏了3p半芯态、或者rc给得太软。
体弹模量也可以作为辅助指标。Ni的实验体弹模量大约在186-190 GPa范围。PBE对这个量的误差普遍比晶格常数大一些，但200 GPa以上就要警惕了。很多情况下体弹模量偏大是“太硬”的赝势表现，这时可以稍微加大rc或者减少投影通道试试。
4.3 截断能收敛曲线
这个检验不值得跳过。用你生成的Ni赝势分别以300、350、400、450、500、600 eV的截断能计算fcc Ni晶胞的总能量（记得固定晶格常数），然后看能量差。
一个合格的PAW_PBE Ni赝势，在400 eV附近应当已经进入能量平台区，500 eV到600 eV之间每原子能量变化应该小于0.01 eV。如果你的赝势在600 eV下还在明显漂移，那多半是rc取得太小，或者径向网格精度不够。
很多新手的误区是“反正PAW硬，截断能设很高就好了”，可高截断能会拖慢计算，又不带来准确度。一张收敛曲线能告诉你的不是“多少eV才能算准”，而是“在多少eV之后增加截断能纯属浪费机时”。
4.4 转移性检验：从金属到表面再到氧化物
真正高质量的PaW_PBE赝势不应该只在一个体系里表现好。Ni催化研究经常横跨fcc金属、表面、纳米颗粒、氧化物界面，你的赝势最好在这些场景里都保持稳定。
转移性检验最实用的做法是跑两个额外体系：
- Ni(111)表面：看表面总能是否明显低于同原子数的随机构型，算一下表面能和文献值对比；
- NiO或Ni与CO的吸附模型：看吸附能是否在合理的1-2 eV范围内波动。
表面上吸附体系的误差如果超过0.1 eV，影响到催化趋势判断时，再去回头调赝势参数。通常这时候你会发现在金属表面计算“看起来不错”的赝势，在吸附体系里会暴露3d投影通道不足的问题。
5. 常见错误、报错排查与实操避坑
5.1 电子构型翻转 / 自洽不收敛
Ni原子计算里最常见的坑就是3d和4s轨道能量接近，自洽过程中电子占据容易在3d和4s之间跳来跳去。这会导致全电子原子计算的“配置翻转”，最终生成赝势里的orbital occupation和你的预期不一致。
出现这种情况时，建议在输入文件里给每个态都填上明确的初始能量，不要用0.0自动猜测。比如3d填0.3、4s填0.1之类的粗糙参考值，至少让程序从接近目标解的轨道出发。另一个办法是把初态配置改成3d9 4s1，这个构型自洽收敛明显更稳。
5.2 投影通道数到底该不该加
有一个容易被人误解的倾向：投影函数越多越准。实际经验是，对Ni来说，s通道一个、p通道一个、d通道两个已经能覆盖绝大部分常规DFT计算。只有在非常精细的电子谱、强关联修正（DFT+U）或者自旋极化明显的体系中，d通道加到3个才能看出微弱好处。
加投影通道的代价不只是计算时间，还可能引入数值近似问题。投影矩阵如果病态，在DFT代码里会导致电子步难以收敛、甚至出现负电子密度。我遇到输出文件里出现“ghost state”或者能带里出现明显伪带的，多数都是投影通道设置过头或者rc与投影函数不匹配导致的。
5.3 磁性计算的初始磁矩
用新生成的Ni赝势算磁性体系时，第一性原理的初始磁矩设置也很关键。对fcc Ni，在DFT输入里设置初始磁矩约0.7 μB/原子，通常能很快收敛到物理态。如果你设置初始磁矩0，可能会掉进顺磁解，磁矩收敛成0，然后你又开始怀疑赝势有问题。其实赝势没问题，是初始磁矩卡在了错误的势阱里。
此外，自旋极化计算里要对s和d通道分开看磁矩贡献。Ni的磁矩主要来自3d通道，4s的磁矩贡献很小，如果生成的赝势在4s上也显示出很强的自旋分裂，那可能是投影通道的s通道设置有问题。
5.4 与官方赝势对比时的预期差
很多人手搓完一套PAW_PBE赝势，立刻拿它跟VASP自带或者pslibrary的Ni赝势对比，发现总能量差了几十甚至上百eV，马上觉得自己生成错了。
这个判断是错的。不同赝势的截止势定义、芯态处理、补偿电荷形式不一样，总能量的绝对零点约定也不同。真正应该对比的不是总能量，而是能量差量——晶格常数、磁矩、吸附能、能带宽度这些相对量。我的习惯是先把两种赝势算出来的物理量放一张表里，如果晶格常数差在0.01 Å以内、磁矩差在0.02 μB以内，这个赝势就是合格的。绝对总能差多少根本不用管。
5.5 后处理转换失败的常见原因
genpaw或转换脚本报错，十有八九不是脚本坏了，而是你安装ATOMPAW时没有把目标代码的支持模块编译进去。也就是说你在configure时没有启用VASP/QE/ABINIT的输出接口，后处理工具自然不认你的目标格式。
处理办法是回到configure阶段，重新检查开关选项。如果确实忘记了，也不用重新生成原子数据，把中间文件存好，重新编译一个有对应接口的ATOMPAW再转换就行。原子数据本身不依赖这些接口，不需要重算。
5.6 一个容易被忽略的硬件相关性
ATOMPAW在部分老版本里对BLAS库的线程数有“怪癖”——你用了多线程OpenBLAS，反而会因为在很小的原子计算中频繁创建线程开销大而变慢，甚至偶发非确定性结果。跑单原子级别计算时，我通常建议把线程数设置为1，比如在运行前设置环境变量：
export OMP_NUM_THREADS=1这样不仅稳定，速度也可能反而更快。等到做大规模晶体计算、真正需要多节点并行时，再把MPI/OpenMP拉满。
几个真正影响Ni赝势成败的细节
最后再聊几个我实际操作中踩过、后来觉得非常值得分享的细节。
第一个是径向网格的“尾端”。很多人看rc小于2.0 bohr，就以为径向网格到2.2 bohr足够，实际上PAW数据集里有一部分形状函数和补偿电荷需要延伸到rc之外一段距离，通常要到2.5-3.0 bohr。截断得太早，会在实空间网格化时产生电荷振荡，表现为总能出现莫名其妙的微小跳跃。
第二个是参考态能量填0.0并不总是安全。虽然ATOMPAW设计了自动猜测机制，但Ni这种3d-4s竞争强烈的元素，我建议至少有一次尝试用手动给定能量。哪怕只是加一个非常粗糙的估计值，自由度变小，自洽路径反而更稳定。
第三个是验证时不要只看一个体系。哪怕你的fcc Ni算得很漂亮，如果打算后续做表面和吸附，最好从一开始就多跑一个表面体系检验。因为表面体系对d投影通道的要求通常更高，等大批量生产计算做完才发现赝势在表面问题上不可靠，那才叫真难受。
生成PAW_PBE赝势这件事，最大的成就感不是“我跑通了”，而是你能清楚说出你这套赝势在物理上做了哪些取舍：参考态是什么、半芯态含不含、rc为什么取这个值、投影通道配到几级。这些决定，恰恰是商业赝势库文档里经常语焉不详的地方。你自己生成一遍，手里才真正有数。
