GPAW无轨道DFT操作

本篇演示GPAW的无轨道DFT（OFDFT）：先生成与所选动能泛函匹配的专用势，再计算原子或周期模型，最后提取能量和密度。

OFDFT把动能写成密度的近似泛函。GPAW这里求解的是与密度平方根有关的集体变量，所以一个计算中的“轨道”可以承载全部电子；这个变量的占据数与普通Kohn–Sham单电子轨道的占据含义不同。

1 下载和运行顺序#

下载OFDFT实操包，保存为gpaw-ofdft.zip。从登录节点提交时，解压后使用sbatch job.slurm；已进入分配好的CPU资源时，依次执行：

python3 -m zipfile -e gpaw-ofdft.zip .
cd gpaw-ofdft
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

任务使用1 MPI进程、1线程。程序先在本算例目录内生成两份专用势，再运行两个条件：λ=2的非周期碳原子，以及λ=1的周期C₈。原生输出、自己生成的XML势和GPW文件都保留在results/；汇总为results/VALIDATION.json，图片在plots/。

2 明确动能泛函和电子数#

本例采用Thomas–Fermi动能与λ倍von Weizsäcker动能的组合，同时使用LDA交换，不加入相关项：

T[n] ≈ T_TF[n] + λ T_W[n]。

在GPAW的这条实现中，Thomas–Fermi部分与交换项通过下面的名称指定：

xcname = '1.0_LDA_K_TF+1.0_LDA_X'

λ还需要传给相应的势生成器和求解器。专用势生成时采用非相对论设置、1.2 Bohr增强球半径、gpernode=150。

本例C的6个电子全部包含在OFDFT势中，冻结芯层数为0；因此一个C原子的电子数是6，C₈是48。它与普通碳PAW设置常见的“4个显式价电子”划分不同，应以当前专用势为准。

3 先生成专用势#

以λ=2为例，下面展示关键操作。文件只写到当前算例自己的子目录：

from pathlib import Path
import os
from gpaw.atom.generator import Generator
from gpaw import setup_paths

folder = Path('setup-lambda2').resolve()
folder.mkdir()
old = Path.cwd()
os.chdir(folder)
try:
    g = Generator('C', xcname=xcname, scalarrel=False,
                  orbital_free=True, tw_coeff=2.0,
                  gpernode=150, txt='generator.txt')
    g.run(rcut=1.2, exx=False, name='lambda_2.0', write_xml=True)
finally:
    os.chdir(old)
setup_paths.insert(0, str(folder))

附件也生成λ=1的对应势。orbital_free=True表示生成OFDFT专用势；tw_coeff与随后计算采用的λ保持一致。程序通过本次进程的setup_paths读取这些文件，不修改集群公共PAW数据集。

4 碳原子：legacy接口与λ缩放#

原子位于10×10×10 Å非周期盒子中心，网格目标间距0.18 Å：

from ase import Atoms
from gpaw import GPAW
from gpaw.eigensolvers import CG
from gpaw.mixer import Mixer

atoms = Atoms('C', positions=[[5, 5, 5]],
              cell=[10, 10, 10], pbc=False)
atoms.calc = GPAW(
    legacy_gpaw=True, mode='fd', h=0.18, xc=xcname,
    setups={'C': 'lambda_2.0'},
    eigensolver=CG(tw_coeff=2.0), mixer=Mixer(0.3, 5, 1),
    convergence={'energy': 1e-6, 'density': 1e-6,
                 'eigenstates': 1e-7}, maxiter=240)
energy = atoms.get_potential_energy()
atoms.calc.write('state.gpw', mode='all')

setups使用刚生成的势名称字符串，便于GPW保存可恢复的标识。保留对应XML文件，并在重启前将它的目录加入setup_paths。

5 周期C₈：新接口小模型#

第二个条件使用8原子的金刚石型立方模型，晶格常数2.8 Å，周期边界，网格目标间距0.14 Å，λ=1：

from ase.build import bulk

atoms = bulk('C', 'diamond', a=2.8, cubic=True)
atoms.calc = GPAW(
    mode='fd', h=0.14, xc=xcname,
    setups={'C': 'lambda_1.0'},
    eigensolver='ppcg', mixer=Mixer(0.1, 5),
    convergence={'energy': 1e-6, 'density': 1e-6,
                 'eigenstates': 1e-7}, maxiter=200)
energy = atoms.get_potential_energy()

这两个条件分别展示当前版本的legacy CG缩放路径与新FD接口路径。C₈的晶格常数是本例输入，结构保持固定；这里没有预测碳材料的平衡晶格常数或实验性质。

6 集群实际结果#

原子能量接近零，是因为GPAW输出使用专用势生成时的原子参考。参考随动能泛函、λ和生成势而定。若要计算材料内聚能，应为材料与孤立原子使用同一套方法和相同λ，再检查数值收敛；不能直接将本表两行代入内聚能公式。

GPAW 26.7.0对应的原生小例给出的近似参考值分别为+0.016和−224.206 eV，本例与它们相符。这是操作与数值流程的对照，不代表这个近似动能泛函能准确描述真实碳材料。

7 密度数据、导出与画图#

图来自本次实际计算的全电子密度。完整数组、晶胞、原子位置、占据和体素体积保存在每个条件的fields.npz中。

source env.sh
"$YS_PY" - <<'PYCODE'
from pathlib import Path
import csv, json
v = json.loads(Path('results/VALIDATION.json').read_text())
with open('ofdft.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['condition', 'lambda', 'energy_eV',
                'electrons', 'collective_occupation'])
    for r in v['results']:
        w.writerow([r['name'], r['lambda_coeff'], r['relative_energy_eV'],
                    r['AE_electrons'], r['occupations'][0]])
        print(r['name'], 'occupation =', r['occupations'])
PYCODE

只重画附件原始数据时：

bash setup.sh
bash plot.sh --reference
ls reference-plots

图片在reference-plots/，不会启动新的势生成或三维自洽计算。下载包提供脚本和自己的提取数据；专用势、GPW检查点由用户运行时生成。

8 环境与使用范围#

本例使用GPAW 26.7.0、Python 3.11.16，CPU构建采用GCC 13.3、HPC-X Open MPI 5.0.10rc2、MKL 2023.2、Libxc 7.0.0。setup.sh加载gpaw/26.7.0-8173m。

OFDFT的结果依赖动能密度泛函及专用势，应用到自己的体系时需要重新检查。需要普通Kohn–Sham轨道、能带、DOS或相应激发分析时，应采用匹配的电子结构方法；这里的单个集体密度变量不构成这些性质的数据。

方法、势生成和λ设置：GPAW OFDFT官方说明。