CP2K 二维晶胞与真空约束

二维层状材料常需要优化面内晶格和原子位置，同时保留真空方向。本篇以三原子单层MoS₂演示固定c的CELL_OPT，明确盒长、晶胞基矢与压力分量各代表什么。

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

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

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

2 先检查六角晶胞的写法#

本例CP2K的HEXAGONAL晶胞采用a、b夹角60°：A=(3.10,0,0)、B=(1.55,2.68467875,0)、C=(0,0,21.14) Å。若从120°六角表示导入，可以用 B_new=B_old+A 改换等价基矢，保持笛卡尔原子坐标不变；它不改变二维Bravais晶格。若使用分数坐标，则需一起正确转换坐标，不能只改B而照抄原分数坐标。

CP2K对CELL_OPT要求晶胞按下三角约定准备，本例A沿x、B在xy面、C沿z。把基矢约定先写清楚，再判断KEEP_ANGLES是否保持了想要的角度。

3 完整输入：固定c，优化面内#

本例用PBE、Mo GTH-PBE-q14和S GTH-PBE-q6、DZVP-MOLOPT-SR-GTH、600/60 Ry、均匀3×3×1 k点。三个原子都允许移动，面内a/b联动保持六角形状。

完整 input.inp：mos2-fixedc-cell-k3#

&GLOBAL
 PROJECT geometry
 RUN_TYPE CELL_OPT
 PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
 METHOD QS
 STRESS_TENSOR ANALYTICAL
 &DFT
  &KPOINTS
   SCHEME MONKHORST-PACK 3 3 1
  &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 8
   &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 3.100000000000 0.000000000000 0.000000000000
   B 1.550000000000 2.684678751732 0.000000000000
   C 0.000000000000 -0.000000000000 21.140000000000
   PERIODIC XYZ
   SYMMETRY HEXAGONAL
  &END CELL
  &COORD
   Mo 0.000000000000 0.000000000000 10.570000000000
   S 1.550000000000 0.894892917244 12.140000000000
   S 1.550000000000 0.894892917244 9.000000000000
  &END COORD
  &KIND Mo
   BASIS_SET DZVP-MOLOPT-SR-GTH
   POTENTIAL GTH-PBE-q14
  &END KIND
  &KIND S
   BASIS_SET DZVP-MOLOPT-SR-GTH
   POTENTIAL GTH-PBE-q6
  &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.
  CONSTRAINT Z
 &END CELL_OPT
&END MOTION

CONSTRAINT Z固定c方向晶胞分量，KEEP_ANGLES TRUE保持角度，KEEP_SYMMETRY TRUE保持指定六角晶胞形状。PERIODIC XYZ仍是带真空的三维周期模型，并不是自动启用严格二维孤立边界。

c=21.14 Å是整个盒长，包括MoS₂层厚，不是纯真空厚度。这里固定盒长来避免优化器把真空当成普通体积自由度；正式性质还需做真空收敛。

4 检查允许的应力分量#

grep -E 'GEOMETRY OPTIMIZATION COMPLETED|Pressure deviation.*constraint' work/cases/mos2-fixedc-cell-k3/output.log
grep -A 5 'STRESS| Analytical stress tensor' work/cases/mos2-fixedc-cell-k3/final-static/output.log
cat work/analysis/RESULTS.json

约束c时，应关注允许变化的面内应力；zz应力不一定为零，不能用它要求继续缩放真空。三维应力的数值还依赖含真空的盒长；需要二维面内张力时，应按一致的盒长和单位转换，而不是直接比较不同真空盒的bar数值。

提取结构和画图#

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点以及允许的自由度。这个三原子例子和当前采样用于流程参考，不代表已经找到所有层型、吸附覆盖度或应变条件的最低态。

相关：块体晶胞与压力、原子与固定原子。方向约束与对称性定义见CP2K2026.2 CELL_OPT手册。