GPAW连续介质溶剂化实操
用固定水分子比较真空、介电屏蔽与表面项，导出真实空腔、介电网格和能量分解。
阅读本篇目录
连续介质方法用空间变化的介电环境描述溶剂，不需要在模型中放入大量显式水分子。本例把一个固定水分子放在隐式水环境中，按“真空→介电常数为1的对照→水介电屏蔽→加入空腔表面项”的顺序操作。
1 运行与查看结果#
下载连续介质溶剂化实操包，保存为gpaw-solvation.zip。新SCF使用CPU计算资源；从登录节点提交可使用sbatch job.slurm。在已分配的CPU资源中：
python3 -m zipfile -e gpaw-solvation.zip .
cd gpaw-continuum-solvation
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh
例程为1 MPI进程、4线程。结果在results/VALIDATION.json，实际空腔与介电网格在results/*-FIELDS.npz。setup.sh会自动准备归档绘图所需的原始网格；如果只想重画本文原图，随后运行：
bash plot.sh --reference
原图在reference-plots/，本次计算图在plots/。
2 固定几何与溶剂参数#
实际水分子由ASE构造，O–H距离约0.9686 Å；使用FD网格间距0.25 Å、PBE和约4 Å的分子外空间。所有条件保持同一几何，不进行结构优化。
atoms = molecule('H2O')
adjust_cell(atoms, 4.0, 0.25)
kwargs = get_HW14_water_kwargs()
kwargs['cavity'].volume_calculator = KB51Volume(
compressibility=0, temperature=298.15
)
calc = SolvationGPAW(
legacy_gpaw=True, mode='fd', xc='PBE', h=0.25,
convergence={'energy':1e-6,'density':1e-5,'eigenstates':1e-8},
maxiter=200, **kwargs
)
原生get_HW14_water_kwargs()提供本例的水模型参数：298.15 K、排斥势强度0.180 eV、H半径1.09 Å、体相介电常数78.36和18.4 mN/m的表面张力参数。
这里变量名epsinf指离溶质很远处的体相极限，不是通常光学语境中水的高频介电常数。空腔函数g接近0时属于溶质区域，局部介电常数接近1；g接近1时进入隐式溶剂，线性模型给出epsilon = 1 + (78.36−1)*g。
3 四个条件怎么设置#
| 条件 | 操作 |
|---|---|
| 真空PBE | 普通GPAW、相同FD网格与几何 |
| ε=1，无表面项 | LinearDielectric(epsinf=1)，interactions=[] |
| ε=78.36，无表面项 | 水介电常数，interactions=[] |
| ε=78.36，含表面项 | 完整HW14水参数 |
完整输入集中在run.py。本例的空腔只依赖几何；表面项通过空腔面积加入，因此适合直接核对能量分解。电子收敛设置不因加入表面项而改变。
4 本集群实算数据#
| 条件 | 总能/eV | 相对真空的变化/eV |
|---|---|---|
| 真空PBE | −14.885494 | 0 |
| ε=1，无表面项 | −14.885569 | −0.000075 |
| ε=78.36，无表面项 | −15.265096 | −0.379603 |
| ε=78.36，含表面项 | −15.207545 | −0.322051 |
ε=1对照与真空只相差约0.075 meV。本例空腔面积为50.11263 Ų，排除体积为33.45641 ų，表面项约为+0.057551 eV；程序独立积分实际网格，核对面积、体积、介电场及表面张力×面积。
get_electrostatic_energy()在该接口中包含除附加短程相互作用以外的DFT能量部分；不要把它误认为单独的库仑能。表中的ΔE是这一固定构型连续介质模型的能量差，没有加入振动、热熵、标准态或构型采样修正，不直接等同于完整实验溶剂化自由能。
5 导出能量和网格#
source env.sh
"$YS_PY" - <<'PYCODE'
import json, csv, numpy as np
from pathlib import Path
f = json.loads(Path('results/VALIDATION.json').read_text())
e0 = f['results'][0]['energy_eV']
with open('water-solvation.csv','w',newline='') as out:
w=csv.writer(out);w.writerow(['condition','energy_eV','delta_eV'])
for r in f['results']:
w.writerow([r['name'],r['energy_eV'],r['energy_eV']-e0])
print(r['name'],r['energy_eV']-e0)
a=np.load('results/HW14eps78_surface-FIELDS.npz')
print('cavity shape:',a['cavity_g'].shape)
print('epsilon range:',a['epsilon'].min(),a['epsilon'].max())
np.save('water-cavity.npy',a['cavity_g'])
np.save('water-dielectric.npy',a['epsilon'])
PYCODE
更换溶质时，要检查原子半径、空腔模型、网格与盒大小；不同溶剂也不能只沿用水的介电常数和表面参数。程序为GPAW 26.7.0 CPU构建，GCC 13.3、HPC-X/OpenMPI 5、MKL 2023.2、Libxc 7.0.0，PAW数据26.8.0。模型说明见GPAW连续介质教程。电子响应的另一个小例见有限电场与极化率。