GPAW分子计算的泊松边界处理

孤立分子使用非周期边界时，盒子边界仍会影响静电势。具有偶极矩的体系，可以尝试多极矩修正或单独扩大泊松方程的网格。本篇保持同一个水分子的几何和电子结构参数，实际比较六种条件，练习怎样设置、检查和导出结果。

这里的较大盒子是有限尺寸对照，并非已知的孤立极限。自己的体系仍应检查盒子、网格、基组以及所需物性的精度。

1 下载、提交和查看#

下载水分子泊松边界算例，保存为gpaw-poisson.zip。从登录节点运行时，解压后用sbatch job.slurm提交；已进入分配好的CPU资源时，可以依次执行：

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

程序使用1个MPI进程、1个线程，完成六组小体系自洽计算。结果在results/VALIDATION.json，原生输出和三维场在各条件子目录，图片在plots/。同样输入重复运行时会复用已完成结果；改变输入请在新目录运行。

2 先明确哪些网格被改变#

电子结构计算使用LCAO的dzp基组、PBE、非周期边界。原子位置固定，不做结构优化。小盒子为7.2×7.2×8.4 Å；较大盒子为12×12×13.2 Å。两个盒子的粗网格间距都为0.30 Å。

扩展泊松网格没有把LCAO计算盒子里的原子复制成超胞，也没有扩大电子波函数的计算盒子。它单独处理静电势边界。小盒子的原始细网格是48×48×56，实际密度和输出势场保存在这个网格上。

3 水分子输入和普通边界#

下面给出小盒子结构，O放在盒中心：

import numpy as np
from ase import Atoms
from gpaw import GPAW, FermiDirac, Mixer
from gpaw.poisson import PoissonSolver

cell = np.array([7.2, 7.2, 8.4])
positions = np.array([[0, 0, 0],
                      [0.7586, 0, 0.5043],
                      [-0.7586, 0, 0.5043]]) + cell / 2
atoms = Atoms('OH2', positions=positions, cell=cell, pbc=False)

calc = GPAW(
    legacy_gpaw=True, mode='lcao', basis='dzp', xc='PBE',
    gpts=(24, 24, 28), nbands=6, spinpol=False,
    symmetry='off', occupations=FermiDirac(0),
    poissonsolver=PoissonSolver('fd', eps=1e-12),
    mixer=Mixer(0.2, 5, 50),
    convergence={'energy': 1e-6, 'density': 1e-6,
                 'eigenstates': 1e-7},
    maxiter=200, txt='water.txt')
atoms.calc = calc
energy = atoms.get_potential_energy()
forces = atoms.get_forces()
dipole = calc.get_dipole_moment()  # e·Å

gpts明确指定网格点数，使两个盒子的网格间距保持一致。eps控制泊松方程求解容差，不等同于电子自洽的能量或密度收敛标准。

4 改用多极矩修正#

把上面的poissonsolver替换为：

from gpaw.poisson_moment import MomentCorrectionPoissonSolver

solver = MomentCorrectionPoissonSolver(
    poissonsolver=PoissonSolver('fd', eps=1e-12),
    moment_corrections=4)
# GPAW(..., poissonsolver=solver, ...)

4表示前4个多极矩，即单极和三个方向的偶极；它不表示4次自洽迭代。改成9还包括五个d型分量。本例两种条件均实际计算。

多极矩数越多，结果不一定在有限网格上单调接近某个有限盒子参考值。应比较目标物性的实际变化，而不是只看这个数字大小。

5 单独扩展泊松网格#

直接在128³扩展细网格求解：

from gpaw.poisson_extravacuum import ExtraVacuumPoissonSolver

solver = ExtraVacuumPoissonSolver(
    gpts=(128, 128, 128),
    poissonsolver_large=PoissonSolver('fd', eps=1e-12))

上面默认coarses=0，不传poissonsolver_small。要先粗化一档再求解，则使用：

solver = ExtraVacuumPoissonSolver(
    gpts=(128, 128, 128),
    poissonsolver_large=PoissonSolver('fd', eps=1e-12),
    coarses=1,
    poissonsolver_small=PoissonSolver('fd', eps=1e-12))

此时原生输出会显示64³的大网格求解，这是粗化后的实际网格。不要把它误读为输入中的128³没有生效。扩展网格必须容纳原细网格；更多粗化也需要满足网格与插值条件，不能只为省内存任意增大coarses。

6 集群实际结果#

下表能量差均相对于“大盒普通边界”，使用相同结构、泛函、基组和网格间距：

本例mom4的能量更接近较大盒子，mom9没有继续单调靠近。扩展网格直接求解与粗化一档的能量差约0.00000267 eV；这个数值属于本小模型和网格设置，不能直接推广到其他体系。

各条件保留了实际原子力。这里的水结构没有优化，所以不能把输出力大小当成已经满足结构优化的证据；有需要时应选择边界处理后再优化，并检查最终受力。

7 电子密度和静电势#

左图是实际水分子全电子总密度，中、右图分别是小盒mom4和extra128相对于小盒普通边界的静电势变化，单位eV。三个图使用相同结构与截面；势场差来自完整同网格数组。密度色标下限以下统一显示为最低颜色，原始数组没有改写。

提取接口如下：

rho = calc.get_all_electron_density(gridrefinement=2)
potential = calc.get_electrostatic_potential()

本例全电子密度积分为10；PAW价电子占据和为8，因为O还包含冻结芯层电子。比较电子数时，需要区分全电子密度积分与价电子占据。

8 导出表格并重画图片#

source env.sh
"$YS_PY" - <<'PYCODE'
from pathlib import Path
import json, csv, numpy as np
r = Path('results')
v = json.loads((r / 'VALIDATION.json').read_text())
with open('water-poisson.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['condition', 'energy_eV', 'dipole_z_e_A', 'AE_electrons'])
    for row in v['results']:
        w.writerow([row['name'], row['energy_eV'],
                    row['dipole_e_A'][2], row['AE_electrons']])
d = np.load(r / 'small-mom4/fields.npz')
print('Density shape:', d['AE_total_density_e_A3'].shape)
print('Potential shape:', d['electrostatic_potential_eV'].shape)
PYCODE

只重画附件自带参考结果：

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

原始参考数组位于reference/，自己的计算数组位于results/；每组fields.npz保存密度、势、原子位置、晶胞、力、偶极矩和占据。附件不包含PAW数据集、GPW检查点或程序二进制。

9 环境与选择建议#

使用GPAW 26.7.0 legacy LCAO、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。

操作时先保留一个普通边界与较大盒子的对照，再比较多极矩或扩展网格；同时看能量、力和所需响应量。周期体系、表面偶极修正和带电体系还涉及各自边界约定，不能直接照搬本非周期中性分子设置。

接口参考：GPAW官方高级泊松求解器说明。