ASE调用CP2K：水分子与网格检查
真实调用计算程序，保存输入输出并用位移能量差分检查力的单位和符号，附终端命令、结构图与完整程序。
阅读本篇目录
ASE用统一的Atoms对象保存结构，把能量和力的计算交给calculator。本篇使用真实程序计算三个邻近构型，再把返回的力与能量差分相互核对。下载本组输入、实际结果及绘图程序，解压进入资料包目录。
1. 准备环境和运行#
本例ASE3.29.0调用CP2K2026.2的CPU构建，通过CP2K shell模式交互，不是只写出一个input文件。水分子在12 Å盒中，三个方向均不周期，采用PBE、DZVP-MOLOPT-SR-GTH和GTH-PBE。
module load ase/3.29.0 cp2k/2026.2-8173m
bash run.sh cp2k
脚本使用2个MPI进程×2个OpenMP线程，并设置足够的栈空间。它调用cp2k.psmp -s，在上下文退出时关闭shell。label="water"保持简短，并切到各案例目录；不要把很长的绝对路径作为CP2K项目名。
原始坐标/Å：
O 6.000 6.000 6.000
H 6.000 6.757 6.587
H 6.000 5.243 6.587
将第2原子沿y先偏移0.03 Å，再取−0.005、0、+0.005 Å。基础网格CUTOFF=400 Ry、REL_CUTOFF=50 Ry，SCF=1E-7。ASE的cutoff参数单位是eV，所以代码写400*Rydberg，不能写成裸的400后声称是400 Ry。
2. 本例完整计算程序#
run.sh负责加载环境和复制输入，实际计算逻辑如下。ASE返回能量单位eV、力单位eV/Å，应力通常为eV/ų。
展开完整Python程序
from pathlib import Path
import json,os,numpy as np
from ase.io import write
from ase.build import bulk
from ase import Atoms
B=Path(__file__).resolve().parent
rows=[]
from ase.calculators.cp2k import CP2K
from ase.units import Rydberg
extra="""&FORCE_EVAL
&DFT
&MGRID
REL_CUTOFF 50
&END MGRID
&SCF
EPS_SCF 1.0E-7
SCF_GUESS ATOMIC
&OT
MINIMIZER DIIS
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&END SCF
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&END DFT
&END FORCE_EVAL"""
for label,delta in [('minus',-.005),('center',0),('plus',.005)]:
d=B/label;d.mkdir();a=Atoms('OHH',positions=[[6,6,6],[6,6.757,6.587],[6,5.243,6.587]],cell=[12]*3,pbc=False);a.positions[1,1]+=.03+delta
os.chdir(d)
with CP2K(label='water',command='mpirun -np 2 cp2k.psmp -s',basis_set='DZVP-MOLOPT-SR-GTH',basis_set_file='BASIS_MOLOPT',potential_file='GTH_POTENTIALS',pseudo_potential='GTH-PBE',xc='PBE',cutoff=400*Rydberg,max_scf=100,poisson_solver=None,stress_tensor=False,inp=extra,print_level='MEDIUM') as calc:
a.calc=calc;e=a.get_potential_energy();f=a.get_forces();write(d/'result.xyz',a);rows.append(dict(label=label,delta_A=delta,energy_eV=e,forces_eV_A=f.tolist()))
(B/'RESULTS.json').write_text(json.dumps(rows,indent=2)+'\n')
生成三个目录：minus、center、plus。程序保存RESULTS.json与每个构型的result.xyz；DFT案例还保留引擎原始输入输出。要继续结构优化时，可将已配置的calculator交给ASE优化器，但本例没有把三个单点误称为优化轨迹。
3. 实际结果和独立检查#
| 位移增量/Å | 总能量/eV |
|---|---|
| -0.005 | -468.569708638 |
| +0.000 | -468.567411173 |
| +0.005 | -468.564297909 |
用中心差分检查第2原子的y方向力：
F ≈ −[E(+δ) − E(−δ)] / (2δ)，δ = 0.005 Å
程序返回：-0.541643009 eV/Å
能量差分：-0.541072912 eV/Å
绝对差异：0.00057 eV/Å
运行python check_results.py reference检查附件参考；复算后执行python check_results.py work检查自己的输出。程序分别核对原始输出和ASE结果；有限差分有截断误差，SCF和积分网格也有误差，不能要求两者逐位相同。相互吻合说明本例接口的主要单位与符号一致，不能证明材料精度已经收敛。
4. 为什么还要看总力与网格#
400 Ry中心构型的总力矢量模约0.00623 eV/Å；孤立分子的整体平移不应有真实恢复力，这个残差提示有限积分网格的影响。提高网格进行一个中心点对照：
bash run.sh cp2k-grid800
| 网格CUTOFF/REL_CUTOFF，Ry | 总力矢量模/eV·Å⁻¹ |
|---|---|
| 400/50 | 0.006232 |
| 800/80 | 0.002171 |
提高网格后残差降低约2.87倍，但仍未达到0.001 eV/Å。中心构型总能量改变0.005951 eV，说明不能忽略数值设置。
本例没有证明800 Ry已充分收敛。正式计算应根据目标精度继续检查CUTOFF、REL_CUTOFF、SCF以及盒尺寸；不要把三个能量差分点当成网格收敛测试。CP2K原始能量是Hartree，接口返回eV，本例按CP2K日志打印的换算常数27.2113838565563 eV/Hartree独立核对；shell已返回eV，不要再乘一次。不同程序内置的物理常数版本可能略有差异，逐位比较时应使用对应程序的常数。
参考：ASE CP2K接口。其他联用：VASP、ABACUS、CP2K、NEP。