鸢算科技 yuansuan Technology
用户门户
CP2K专题

CP2K温度加速Monte Carlo：辅助变量与接受率

使用真实水二聚体TAMC，区分物理与辅助温度，比较原子步长并读取SCF、接受率和广义力。

更新于 2026-11-12约 11 分钟阅读按任务与环境选择适用配置
阅读本篇目录

本系列用鸢算科技集群实际完成的三个CP2K入口，说明传统Monte Carlo、Tree Monte Carlo和温度加速Monte Carlo的输入、提交与结果读取。它们解决的问题和计数方式不同，不能只改RUN_TYPE名字就互换。

下载完整输入、实际输出与绘图工具，保存为example.zip：

python3 -m zipfile -e example.zip .
cd cp2k-special-sampling-examples
bash results.sh
bash read.sh
bash plot.sh

以上只读已经完成的结果，不启动新计算。程序为cp2k/2026.2-8173m，Intel oneAPI 2024.2/ifx、Intel MPI 2021.13、MKL构建。原始方法输入改编自固定CP2K 2026.2官方示例，附GPL许可；三种例子均为小规模操作参考，不据此报告相平衡、扩散时间或收敛自由能。

要重算整套小例子，在新的解压目录提交：

sbatch run.slurm
squeue -u "$USER"
# 完成后
bash read.sh work analysis-new
bash plot.sh analysis-new figures-new work

CPU入口按集群规则申请整节点56核，实际按每个步骤分别用1、4或2个MPI进程，OpenMP为1；这不是56核性能测试。work/已经存在时脚本会保留它，避免覆盖结果。

1. 物理温度与辅助温度分开#

Temperature Accelerated Monte Carlo引入辅助集体变量，以耦合势连接原子坐标和辅助变量。本例是六原子水二聚体，HCTH120/DZVP-GTH、300 Ry、电子SCF阈值1E−6；物理系统温度25 K，辅助变量温度3000 K。这是沿用官方小例的算法演示，不表示水真的同时具有两种物理温度，也不是通用推荐参数。

集体变量用原版输入的HBP定义，ATOMS 4 5 1中的编号对应水二聚体文件；它不是一个可以直接标成Å的普通键长。辅助坐标z、原子构型给出的Θ和广义力要按该定义解释。

cat inputs/water-dimer-TAMC-small-step/dimer.xyz
cat inputs/water-dimer-TAMC-small-step/input.inp
完整TAMC输入：展开查看DFT、HBP、辅助变量和MC/MD控制 
#CPQA INCLUDE dimer.xyz
#CPQA INCLUDE topology_atoms.psf
!!cp2k!!
!! inp file by Alin M Elena
!! alinm.elena@gmail.com
!! http://alin.elenaworld.net
#global section sets what kind of calculations you want to do an general
! settings
@SET pfreq 100 ! the frequency for printing .cube properties and restart files
@SET name dimer
@SET period T
@SET method GPW ! use GAPW/ GPW for dft
@SET zk 10.0  ! coupling constant for collective variables and atoms
@SET zmass 100.0 ! mass for collective variables
@SET zgamma 0.001 ! friction coefficient for lagrange equations for cvs
@SET Tz  3000.0 ! temperature
@SET Nz 1   ! number of steps to integrate z
@SET dtz 0.5 ! time step to integrate the z
@SET dt 0.1 ! time steps to integrate atoms
@SET N 1 ! the length of the mini-trajectory
@SET M 2 ! number of cycles/steps for the HMC part
@SET Temp 25 ! physical system temperature
&GLOBAL
  ECHO_INPUT yes
  PREFERRED_FFT_LIBRARY FFTSG
  !   Run_type MD
  !  Run_type GEOMETRY_OPTIMIZATION
  PRINT_LEVEL low
  PROJECT ${name}
  !   Run_type ENERGY_Force
  RUN_TYPE TAMC
  WALLTIME 3600
  &PRINT
    FILENAME =./${name}
  &END PRINT
&END GLOBAL

#defines how forces are evaluated
&FORCE_EVAL
  METHOD QS  ! what is the method used for computing the energy
  STRESS_TENSOR analytical
  &DFT  ! dft settings
    BASIS_SET_FILE_NAME BASIS_SET
    CHARGE 0
    POTENTIAL_FILE_NAME POTENTIAL
    &LOCALIZE True
      EPS_LOCALIZATION 1.0e-4
      EPS_OCCUPATION 1.0e-6
      MIN_OR_MAX SpreadMin
      RESTART false
      STATES Occupied
      USE_HISTORY true
      &PRINT
        &MOLECULAR_DIPOLES Off
          ADD_LAST numeric
          PERIODIC ${period}
          REFERENCE zero
          REFERENCE_POINT 0.0 0.0 0.0
          &EACH
            MD 10
          &END EACH
        &END MOLECULAR_DIPOLES
        &MOLECULAR_STATES off
          ADD_LAST numeric
          &EACH
            MD 10
          &END EACH
        &END MOLECULAR_STATES
        &TOTAL_DIPOLE On
          ADD_LAST numeric
          PERIODIC ${period}
          REFERENCE zero
          REFERENCE_POINT 0.0 0.0 0.0
          &EACH
            MD 10
          &END EACH
        &END TOTAL_DIPOLE
        &WANNIER_CENTERS
          ADD_LAST numeric
          FORMAT xyz
          IONS+CENTERS true
          &EACH
            MD 1
          &END EACH
        &END WANNIER_CENTERS
        &WANNIER_STATES on
          ADD_LAST numeric
          &EACH
            MD 1
          &END EACH
        &END WANNIER_STATES
      &END PRINT
    &END LOCALIZE
    &MGRID
      CUTOFF 300
    &END MGRID
    &POISSON
      PERIODIC None
      POISSON_SOLVER MT
      &MT
      &END MT
    &END POISSON
    &PRINT
      &EFIELD_CUBE Medium
        ADD_LAST Numeric
        &EACH
          MD ${pfreq}
        &END EACH
      &END EFIELD_CUBE
      &E_DENSITY_CUBE Medium
        ADD_LAST Numeric
        &EACH
          MD ${pfreq}
        &END EACH
      &END E_DENSITY_CUBE
      &MO_CUBES Medium
        ADD_LAST Numeric
        &EACH
          MD 1
        &END EACH
      &END MO_CUBES
      &TOT_DENSITY_CUBE Medium
        ADD_LAST Numeric
        &EACH
          MD ${pfreq}
        &END EACH
      &END TOT_DENSITY_CUBE
      &V_HARTREE_CUBE Medium
        ADD_LAST Numeric
        &EACH
          MD ${pfreq}
        &END EACH
      &END V_HARTREE_CUBE
    &END PRINT
    &QS
      EPS_DEFAULT 1.0e-12
      EXTRAPOLATION ASPC
      EXTRAPOLATION_ORDER 3
      METHOD ${method}
      &DFTB
        DISPERSION T
        DO_EWALD T
        ORTHOGONAL_BASIS F
        SELF_CONSISTENT T
        &PARAMETER
          PARAM_FILE_NAME scc_parameter
          PARAM_FILE_PATH params_scc
          UFF_FORCE_FIELD uff_table
        &END PARAMETER
      &END DFTB
    &END QS
    &SCF  ! the scf settings...
      EPS_SCF 1.e-6
      MAX_SCF 100
      SCF_GUESS Atomic
      &MIXING
        METHOD Pulay_Mixing
        NBUFFER 5
        NMIXING 1
        PULAY_ALPHA 0.3
        PULAY_BETA 0.5
      &END MIXING
      &OT T  ! orbital transformation settings
        MINIMIZER CG
        PRECONDITIONER Full_single_inverse
      &END OT
      &OUTER_SCF On
        EPS_SCF 1.0E-6
        MAX_SCF 10
      &END OUTER_SCF
      &PRINT
        &RESTART
          ADD_LAST Numeric
          &EACH
            GEO_OPT 1
            MD ${pfreq}
            QS_SCF 0
          &END EACH
        &END RESTART
      &END PRINT
    &END SCF
    &XC  ! exchange and correlation functionals
      &XC_FUNCTIONAL HCTH120
        &HCTH T
          PARAMETER_SET 120
        &END HCTH
      &END XC_FUNCTIONAL
    &END XC
  &END DFT
  &SUBSYS
    &CELL
      ABC [angstrom] 6.0 6.0 6.0
      ALPHA_BETA_GAMMA 90.0 90.0 90.0
      PERIODIC xyz
    &END CELL
    &COLVAR
      &HBP
        ATOMS 4 5 1
        NPOINTS 1
      &END HBP
    &END COLVAR
    &KIND O
      BASIS_SET DZVP-GTH-PADE
      POTENTIAL GTH-HCTH120-q6
    &END KIND
    &KIND H
      BASIS_SET DZVP-GTH-PADE
      POTENTIAL GTH-HCTH120-q1
    &END KIND
    &TOPOLOGY
      CONNECTIVITY MOL_SET
      COORD_FILE_FORMAT xyz
      COORD_FILE_NAME dimer.xyz
      &CENTER_COORDINATES off
      &END CENTER_COORDINATES
      &MOL_SET
        &MOLECULE
          CONN_FILE_NAME topology_atoms.psf
          NMOL 1
        &END MOLECULE
      &END MOL_SET
    &END TOPOLOGY
  &END SUBSYS
&END FORCE_EVAL

&MOTION
  &FREE_ENERGY Metadyn
    &METADYN
      LAGRANGE true
      LANGEVIN true
      TAMCSTEPS ${Nz}
      TEMPERATURE ${Tz}
      TEMP_TOL 2000.0
      TIMESTEP ${dtz}
      &EXT_LAGRANGE_SS0
        0.44025202
      &END EXT_LAGRANGE_SS0
      &METAVAR
        COLVAR 1
        GAMMA ${zgamma}
        LAMBDA ${zk}
        MASS ${zmass}
        SCALE 1.0
      &END METAVAR
      &PRINT
        &COLVAR
          &EACH
            MD 1
          &END EACH
        &END COLVAR
      &END PRINT
    &END METADYN
  &END FREE_ENERGY
  &GEO_OPT
    MAX_DR 3.0e-3
    MAX_FORCE 4.5e-4
    MAX_ITER 1000
    OPTIMIZER LBFGS
  &END GEO_OPT
  &MC
    ENSEMBLE Traditional
    IPRINT 1
    NSTEP ${M}
    RANDOMTOSKIP 0
  &END MC
  &MD
    ENSEMBLE NVE
    STEPS ${N}
    TEMPERATURE ${Temp}
    TIMESTEP ${dt}
    TIME_START_VAL 0.0
    &THERMOSTAT
      &NOSE
        LENGTH 3
        MTS 2
        TIMECON 1000.0
        YOSHIDA 3
      &END NOSE
    &END THERMOSTAT
  &END MD
  &PRINT
    &TRAJECTORY
      FILENAME trajectory.xyz
      FORMAT xyz
      &EACH
        MD 10
      &END EACH
    &END TRAJECTORY
    &VELOCITIES
      FILENAME velocity.xyz
      FORMAT xyz
      &EACH
        MD 10
      &END EACH
    &END VELOCITIES
  &END PRINT
&END MOTION

关键变量中，dt=0.1是原子小轨迹的时间步，dtz=0.5是辅助变量积分步；M=2为每个阶段的MC循环数，Nz=1只推进一个辅助步。初始平均阶段还要执行两次小轨迹，所以本次共尝试四条小轨迹。不要把四条尝试直接当作四个接受状态或一条常规MD轨迹。

2. 运行并观察电子收敛#

在新的目录保留input.inp、dimer.xyz与topology_atoms.psf，并在已申请的CPU资源中运行：

module purge
module load cp2k/2026.2-8173m
export OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1
mpirun -np 2 cp2k.psmp -i input.inp -o output.log
grep 'SCF run converged' output.log
grep 'HMC| global acceptance ratio' output.log

本次每个TAMC条件各完成9次SCF受力评估。它们是采样内部的电子求解，不是9个已经热化的独立材料样本。

3. 为什么比较两个原子时间步#

附件保留原子时间步1 fs和0.1 fs的两组实际结果，其他模型与温度参数一致。1 fs例四条小轨迹均拒绝；仅把原子步长改为0.1 fs后，初始两条接受，后两条拒绝，累计2/4。

cat analysis-read/CP2K-TAMC-step-comparison.csv
cat analysis-read/CP2K-TAMC-native-forces.csv

实际时间步与接受率对照、辅助广义力

在初始固定z平均阶段，耦合广义力与−10(z−Θ)对应，配套读取器核对实际数值与运行平均。后续辅助变量更新阶段涉及算法中的平均力和状态更新，不把每一列都机械当作同一时刻的瞬时原子力。

这个对照说明需要关注积分误差和接受率；四次尝试太少，2/4不是可靠的长期接受率估计，更不能据此给出自由能曲线或反应速率。

4. 改成自己的研究前#

先定义物理体系、集体变量、耦合强度、辅助质量和两种温度的角色，再检查电子求解、积分步、热化和足够长的统计。需要材料自由能时，要按方法定义进行重构并检查收敛，不能将本例两行辅助变量输出直接画成PMF。

下载时间步对照表 · 下载原生广义力记录。更多参数见CP2K 2026.2 MC输入、TMC输入及辅助变量METADYN设置。

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

查看配图

放大后可滚动查看细节，按 Esc 关闭。

100%查看原图 ↗