CP2K 晶胞优化与压力

想寻找给定外压下的晶胞参数，应使用 RUN_TYPE CELL_OPT。它允许晶胞和原子坐标一起优化；对称性、角度、体积或方向约束决定哪些晶胞自由度可以变化。本篇用Si8立方晶胞，保留立方形状，调整边长。

下载完整输入、真实结构和原始输出、数据与绘图脚本，解压后进入 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 cell

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

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

2 指定立方晶胞与优化目标#

初始a=5.35 Å，PBE/DZVP-MOLOPT-SR-GTH、Si GTH-PBE-q4、600/60 Ry、2×2×2 k网格、对角化与Broyden混合。STRESS_TENSOR ANALYTICAL 请求解析应力。

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

&GLOBAL
 PROJECT geometry
 RUN_TYPE CELL_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.350000000000 0.000000000000 0.000000000000
   B 0.000000000000 5.350000000000 0.000000000000
   C 0.000000000000 0.000000000000 5.350000000000
   PERIODIC XYZ
   SYMMETRY CUBIC
  &END CELL
  &COORD
   Si 0.000000000000 0.000000000000 0.000000000000
   Si 1.337500000000 1.337500000000 1.337500000000
   Si 0.000000000000 2.675000000000 2.675000000000
   Si 1.337500000000 4.012500000000 4.012500000000
   Si 2.675000000000 0.000000000000 2.675000000000
   Si 4.012500000000 1.337500000000 4.012500000000
   Si 2.675000000000 2.675000000000 0.000000000000
   Si 4.012500000000 4.012500000000 1.337500000000
  &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
 &CELL_OPT
  OPTIMIZER BFGS
  MAX_ITER 80
  MAX_FORCE 3.5E-4
  RMS_FORCE 2.0E-4
  MAX_DR 0.003
  RMS_DR 0.0015
  EXTERNAL_PRESSURE 0
  PRESSURE_TOLERANCE 200
  KEEP_ANGLES .TRUE.
  KEEP_SYMMETRY .TRUE.
 &END CELL_OPT
&END MOTION

关键设置：SYMMETRY CUBIC 表明立方晶胞，KEEP_SYMMETRY TRUE 保持立方形状；KEEP_ANGLES TRUE 保持角度，EXTERNAL_PRESSURE 0 以零外压为目标，压力单位bar；容差200 bar=0.02 GPa。这里没有设置KEEP_VOLUME，因此允许体积变化。

CELL_OPT同时放松原子。若保持体积或固定某些方向，要明确增加相应约束，不能把无约束块体优化直接用于真空薄层。

3 读取晶格、力与应力#

grep -E 'GEOMETRY OPTIMIZATION COMPLETED|Internal pressure|Pressure is converged' work/cases/si-cell-opt-k2/output.log
grep -A 5 'STRESS| Analytical stress tensor' work/cases/si-cell-opt-k2/final-static/output.log
cat work/analysis/RESULTS.json

压力由实际应力张量的迹除以三取得。原生日志还会列出应力本征向量，不要把后面的本征向量表当成应力张量。

4 用最终结构做后续性质#

脚本在 final-static 中生成以最终原子与晶胞为输入的 ENERGY_FORCE，重新自洽求能量、力和应力。后续DOS、能带或声子应从确认后的最终结构开始，不要仍使用初始5.35 Å晶胞。

提取结构和画图#

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。

本例用于学习流程，给出的a依赖有限高斯基组、网格、赝势和k采样；需要材料参数时，应继续做这些设置的收敛。只要原子位置变化、晶胞不变的情况看原子与约束优化；薄层看二维晶胞。设置依据CP2K2026.2 CELL_OPT手册。