CP2K 原子位置与固定原子优化

想保持晶格不变，只放松原子时，用 RUN_TYPE GEO_OPT。本篇以Si8为例，从a=5.43 Å金刚石结构出发，给第二个原子一个小位移，再分别演示全部原子可移动和固定第一个原子。

下载完整输入、真实结构和原始输出、数据与绘图脚本，解压后进入 run.sh 所在目录。采用CP2K 2026.2 P100构建：GNU13.3、MPICH5.0.1、CUDA12.9、OpenBLAS0.3.33、ScaLAPACK2.2.3，双精度 cp2k.psmp；每个工作流使用一个MPI进程、四个OpenMP线程和一张P100。分析用ASE3.29.0。

1 提交这一组#

sbatch submit.slurm atoms

这是核时P100模板。套餐用户用 yskit 生成自己的分区脚本，保留脚本里的 bash run.sh atoms；在分配的计算资源中执行，不在登录节点运行科学计算。

inputs 是可修改输入，reference/cases 是本文实际结果，work/cases 是你的新计算。脚本按“优化→取最终结构→静态单点→检查力/应力→转换结构→画图”顺序执行，不只停在程序退出。已有work时请新解压，避免覆盖原结果。

2 准备结构和电子设置#

&COORD中的坐标为Å，本例晶胞固定为5.43 Å立方。使用PBE、DZVP-MOLOPT-SR-GTH、Si的GTH-PBE-q4、600 Ry网格与REL_CUTOFF60 Ry，均匀2×2×2 k点；采用对角化与Broyden混合，EPS_SCF为1E−7。两次优化的晶胞、赝势和采样相同。

完整 input.inp：si-atom-opt-k2#

&GLOBAL
 PROJECT geometry
 RUN_TYPE GEO_OPT
 PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
 METHOD QS
 STRESS_TENSOR ANALYTICAL
 &DFT
  &KPOINTS
   SCHEME MONKHORST-PACK 2 2 2
  &END KPOINTS
  BASIS_SET_FILE_NAME ../DATA/BASIS_MOLOPT
  POTENTIAL_FILE_NAME ../DATA/GTH_POTENTIALS
  &QS
   EPS_DEFAULT 1.0E-12
  &END QS
  &MGRID
   CUTOFF 600
   REL_CUTOFF 60
  &END MGRID
  &SCF
   EPS_SCF 1.0E-7
   MAX_SCF 120
   SCF_GUESS ATOMIC
   ADDED_MOS 4
   &DIAGONALIZATION
    ALGORITHM STANDARD
   &END DIAGONALIZATION
   &MIXING
    METHOD BROYDEN_MIXING
    ALPHA 0.3
    NBROYDEN 8
   &END MIXING
  &END SCF
  &XC
   &XC_FUNCTIONAL PBE
   &END XC_FUNCTIONAL
  &END XC
 &END DFT
 &SUBSYS
  &CELL
   A 5.430000000000 0.000000000000 0.000000000000
   B 0.000000000000 5.430000000000 0.000000000000
   C 0.000000000000 0.000000000000 5.430000000000
   PERIODIC XYZ
   SYMMETRY NONE
  &END CELL
  &COORD
   Si 0.000000000000 0.000000000000 0.000000000000
   Si 1.427500000000 1.317500000000 1.407500000000
   Si 0.000000000000 2.715000000000 2.715000000000
   Si 1.357500000000 4.072500000000 4.072500000000
   Si 2.715000000000 0.000000000000 2.715000000000
   Si 4.072500000000 1.357500000000 4.072500000000
   Si 2.715000000000 2.715000000000 0.000000000000
   Si 4.072500000000 4.072500000000 1.357500000000
  &END COORD
  &KIND Si
   BASIS_SET DZVP-MOLOPT-SR-GTH
   POTENTIAL GTH-PBE-q4
  &END KIND
 &END SUBSYS
 &PRINT
  &FORCES ON
  &END FORCES
  &STRESS_TENSOR ON
  &END STRESS_TENSOR
 &END PRINT
&END FORCE_EVAL
&MOTION
 &GEO_OPT
  OPTIMIZER BFGS
  MAX_ITER 80
  MAX_FORCE 3.5E-4
  RMS_FORCE 2.0E-4
  MAX_DR 0.003
  RMS_DR 0.0015
 &END GEO_OPT
&END MOTION

3 固定原子怎么写#

约束例 si-fixed-opt-k2 使用相同电子和结构设置，只增加：

&MOTION
 &CONSTRAINT
  &FIXED_ATOMS
   LIST 1
  &END FIXED_ATOMS
 &END CONSTRAINT
&END MOTION

实际文件把该段放在已有的同一个MOTION内，不需要在原文件末尾重复建立MOTION。LIST 1 指坐标列表中第一个原子，编号从1开始；不是Si的元素编号。改变结构顺序后必须重新核对固定原子。

几何力目标 MAX_FORCE 3.5E-4 单位Hartree/Bohr，约0.0180 eV/Å；RMS力和位移也参与停止条件。最终静态同时列出最大分量与最大矢量模，便于按自己的标准判断。

4 查看优化和最终单点#

grep -E 'GEOMETRY OPTIMIZATION COMPLETED|Maximum gradient' work/cases/si-*/output.log
grep -E 'ENERGY\||Atomic forces|STRESS\| Analytical' work/cases/si-*/final-static/output.log
cat work/analysis/RESULTS.json

固定原子坐标保持不变，其余原子完成约束下的优化。约束原子的残余受力不应与可移动自由度混在一起判断；改变固定区域会改变优化问题和最终能量。

5 原子力很小，晶胞未必平衡#

本例固定a=5.43 Å之后原子力很小，但最终静态平均压力仍约 2.126 GPa。这是当前有限参数模型下固定晶胞的应力，说明原子位置优化没有自动优化晶胞。不要因为GEO_OPT结束，就认定体积和压力也达到平衡。

提取结构和画图#

module load ase/3.29.0
python analyze.py work
python convert_structures.py work
python plot_results.py work
ls work/converted
ls work/figures

analyze.py读取原生日志，convert_structures.py实际调用 ase convert，将原生最终CIF转换为EXTXYZ，核对元素、坐标和晶胞；图片均从这些结构与数据生成。只重画本文结果时，把命令中的work换成reference。

换成自己的体系时，先确定要固定晶胞还是同时优化晶胞，再检查坐标单位、原子编号、k点和电子收敛。继续看晶胞与压力优化和二维真空约束。