CP2K专题

CP2K力匹配：参考轨迹到势参数拟合

用8原子Ar实际参考轨迹演示OPTIMIZE_INPUT与FORCE_MATCHING，给出51帧数据、11次评估和参数读取。

更新于 2026-11-13约 8 分钟阅读按任务与环境选择适用配置

阅读本篇目录

想理解CP2K如何把参考构型与力用于参数拟合，可以先做一个小的Lennard–Jones例子。本篇先用已知Ar参数产生参考数据，再让OPTIMIZE_INPUT调整两个参数，观察目标函数如何变化。

参考数据来自本例CP2K经典势MD，不是DFT数据；目的在于看清文件传递和命令流程。11次函数评估达到预算后停止，得到的参数不能直接作为新的可靠Ar力场使用。

1. 下载、看结构与数据#

下载CP2K原子、批量与力匹配资料包，保存为example.zip：

python3 -m zipfile -e example.zip .
cd cp2k-atom-farming-forcefit
bash results.sh
bash read.sh
bash plot.sh

这几条命令先读取随包保存的真实计算结果。CSV在data/，图片在figures/，完整输入和原生输出分别在inputs/与reference-results/。

8原子Ar结构和实际参考轨迹

参考设置	本例值
原子数	8 Ar
系综与温度	NPT_I，300 K
步数与步长	100步，2 fs
保存间隔	每2步，含初始帧，共51帧
LJ参数	ε=119.8 K，σ=3.4 Å
截断	12 Å
参考MD随机种子	802101

这是200 fs短流程，不用于平衡性质或统计收敛判断。

2. 第一步：准备参考计算#

cat inputs/force-matching/argon-ref.inp

展开完整参考MD输入

&GLOBAL
  IOLEVEL SILENT
  PROJECT REF
  SEED 802101
  RUN_TYPE MD
&END GLOBAL

&MOTION
  &MD
    ENSEMBLE NPT_I
    STEPS 100
    TEMPERATURE 300
    TIMESTEP 2.0
    &BAROSTAT
    &END BAROSTAT
    &THERMOSTAT
      &NOSE
        LENGTH 3
        MTS 2
        TIMECON 100
        YOSHIDA 3
      &END NOSE
    &END THERMOSTAT
  &END MD
  &PRINT
    &CELL
      &EACH
        MD 2
      &END EACH
    &END CELL
    &FORCES
      &EACH
        MD 2
      &END EACH
    &END FORCES
    &TRAJECTORY
      &EACH
        MD 2
      &END EACH
    &END TRAJECTORY
  &END PRINT
&END MOTION

&FORCE_EVAL
  METHOD FIST
  STRESS_TENSOR ANALYTICAL
  &MM
    &FORCEFIELD
      &CHARGE
        ATOM Ar
        CHARGE 0.0
      &END CHARGE
      &NONBONDED
        &LENNARD-JONES
          ATOMS Ar Ar
          EPSILON 119.8
          RCUT [angstrom] 12.0
          SIGMA [angstrom] 3.4
        &END LENNARD-JONES
      &END NONBONDED
    &END FORCEFIELD
    &POISSON
      &EWALD
        EWALD_TYPE none
      &END EWALD
    &END POISSON
  &END MM
  &SUBSYS
    &CELL
      A 9.6310875301855674E+00 0.0000000000000000E+00 0.0000000000000000E+00
      B 0.0000000000000000E+00 9.6310875301855674E+00 0.0000000000000000E+00
      C 0.0000000000000000E+00 0.0000000000000000E+00 9.6310875301855674E+00
    &END CELL
    &COORD
      Ar   -2.9801516326320010E-01   -1.3601907871663195E+00    5.5171130487752702E+00
      Ar    3.4149510515239569E+00    1.0652576364029100E+00    5.5460851143546490E+00
      Ar   -4.8038262016805466E-01    4.9219660245739698E+00    5.4531759530226092E+00
      Ar    4.3899163615664634E+00    6.3696825498493626E+00    4.7705625410379264E+00
      Ar   -5.2621394328845272E-01    1.1638196685656448E+00   -5.1719754234841919E-01
      Ar    3.9341474087017749E+00   -1.1516822001738987E+00   -3.0662688763800350E-01
      Ar    1.4455908175038696E+00    4.1262301758627968E+00   -4.7434313907940118E-01
      Ar    7.3821811477948271E+00    4.1270919924567533E+00   -7.2659402775337623E-01
    &END COORD
    &TOPOLOGY
    &END TOPOLOGY
  &END SUBSYS
&END FORCE_EVAL

提交脚本中的第一条参考计算是：

cp2k.psmp -i argon-ref.inp > reference.out

得到以下三类文件后才能进入拟合，帧数与步号需要对应：

ls reference-results/force-matching/REF-pos-1.xyz
ls reference-results/force-matching/REF-frc-1.xyz
head reference-results/force-matching/REF-1.cell

位置XYZ的坐标是Å，力XYZ沿用CP2K原生原子单位，晶胞文件是Å。不要把力未经换算替换成eV/Å，再继续沿用同一套拟合输入。

3. 第二步：指定待拟合参数#

argon-fit.inp中的LJ参数写成变量占位符：

EPSILON ${PRM01}
SIGMA [angstrom] ${PRM02}

完整待拟合模板可读：

cat inputs/force-matching/argon-fit.inp

模板中还保留MD段落，但在本例FORCE_MATCHING调用中使用参考帧来评价势，不是在每次评价时重新生产那条100步参考MD。

主拟合输入：

#CPQA DEPENDS argon-ref.inp
#CPQA INCLUDE argon-fit.inp
&GLOBAL
  PROGRAM OPTIMIZE_INPUT
  PROJECT FITIT
  RUN_TYPE NONE
  SEED 1
&END GLOBAL

&OPTIMIZE_INPUT
  ACCURACY 0.0002
  MAX_FUN 11
  METHOD FORCE_MATCHING
  ! will be randomized by 10%
  RANDOMIZE_VARIABLES 10
  STEP_SIZE 0.5
  &FORCE_MATCHING
    GROUP_SIZE 2
    OPTIMIZE_FILE_NAME argon-fit.inp
    REF_CELL_FILE_NAME REF-1.cell
    REF_FORCE_FILE_NAME REF-frc-1.xyz
    REF_TRAJ_FILE_NAME REF-pos-1.xyz
    SHIFT_MM 0.0
    SHIFT_QM 0.0
    &COMPARE_ENERGIES
      &EACH
        POWELL_OPT 5
      &END EACH
    &END COMPARE_ENERGIES
  &END FORCE_MATCHING
  &HISTORY
    &EACH
      POWELL_OPT 5
    &END EACH
  &END HISTORY
  &RESTART
    &EACH
      POWELL_OPT 5
    &END EACH
  &END RESTART
  &VARIABLE
    LABEL PRM01
    VALUE 119.8
  &END VARIABLE
  &VARIABLE
    LABEL PRM02
    VALUE 3.4
  &END VARIABLE
&END OPTIMIZE_INPUT

PROGRAM OPTIMIZE_INPUT与PROGRAM_NAME OPTIMIZE_INPUT在该输入中同义。REF_TRAJ_FILE_NAME、REF_FORCE_FILE_NAME、REF_CELL_FILE_NAME指向刚生成的参考数据；OPTIMIZE_FILE_NAME指向带变量的势模板。GROUP_SIZE 2在本例对应2个MPI进程组成一组。

RANDOMIZE_VARIABLES 10会先扰动初始变量，所以真实起点不是恰好119.8和3.4。MAX_FUN 11限制函数评估次数，不代表执行11轮独立MD，也不保证参数已优化收敛。

4. 第三步：运行并读取拟合#

本例采用集群cp2k/2026.2-8173m，Intel oneAPI 2024.2和Intel MPI构建，CPU路径，不使用GPU。新解压目录中运行：

sbatch run.slurm

脚本依次执行原子、两个FARMING子任务、Ar参考MD和参数拟合。按整节点规则申请56核，最多2个MPI进程、每进程1线程；这是功能教学，不是性能测试。完成后：

bash read.sh work analysis-new
bash plot.sh analysis-new figures-new work

已有结果可直接读取，不必等待CPU队列，也不要把新模拟直接放在登录节点。

拟合阶段的实际命令是：

mpirun -np 2 cp2k.psmp -i driver.inp > fit.out

实际11次拟合评价和历史最优值

grep -A 3 'initial values' reference-results/force-matching/fit-native.log
grep -A 3 'optimal variables found so far' reference-results/force-matching/fit-native.log
cat data/CP2K-native-force-matching-cost.csv

项目	实际结果
扰动后初始ε/K	107.8281362703
扰动后初始σ/Å	3.4379988869
初始目标值	0.000127994501
11次中的最优目标值	0.000083099444
最优ε/K	107.2759501441
最优σ/Å	3.4235368089

目标值降低约35.1%，但输出明确是达到最大评估次数后停止。本例只说明原生拟合链路和读数方法，不能用目标值变小代替参数可靠性验证。

5. 哪个文件保存最终最优参数#

本例HISTORY和RESTART每5次输出，历史文件只包含第5、10次，不包含第11次最优点；重启文件也只代表写出时的状态。本例最终最优变量应从主输出的optimal variables found so far读取，不要把最后一行稀疏历史当成最终最优参数。

自己的研究还需要更丰富的参考构型、独立验证集、适当的目标函数与足够的优化预算；本例由同类LJ模型生成参考，不能代表对真实材料或DFT力的拟合精度。

需要根据你的体系选择资源？联系讨论计算需求。实际可用软件与权限，以当前账户和运行环境为准。