CP2K与i-PI温度副本交换

温度副本交换同时运行不同温度的体系，按Monte Carlo规则交换温度所属的系综，帮助构型访问不同采样条件。CP2K负责计算真实势能和原子力，i-PI负责核坐标推进、恒温与交换。本篇从一个能独立理解的三维谐振子小模型开始操作。

1 下载、运行与查看目录#

下载CP2K与i-PI温度交换完整操作包，保存为cp2k-replica-exchange.zip。普通作业可解压后用sbatch job.slurm提交；下面短命令用于已经通过Slurm分配的CPU计算环境。

CPU分区按56核整节点申请与计费。包内job.slurm已设置单节点、1个Slurm任务、每任务56核；案例脚本在该分配中按各自所需的少量核心运行，每个CP2K进程仍为1线程。

python3 -m zipfile -e cp2k-replica-exchange.zip .
cd cp2k-replica-exchange
bash setup.sh

bash run.sh
bash verify.sh
bash plot.sh

脚本依次启动一个i-PI服务进程和两个CP2K力计算进程，运行4000步，关闭本地连接后调用原生i-pi-remdsort，再提取数据画图。每个进程用一个物理核心，CP2K为1个MPI进程、1线程。结果在results/，本次运行的图在plots/；reference-plots/复现的是操作包内实际参考算例的图片。

CP2K为cp2k/2026.2-8173m，Intel ifx2024.2.1、IntelMPI2021.13、MKL2024.2。操作包携带原版i-PI3.3.1纯Python程序及许可证，setup.sh仅在当前工作目录展开，NumPy使用集群已有环境。

2 本例算什么#

模型为一个质量1 u的粒子，处在U=½K(x²+y²+z²)的三维外部谐振势中，K=0.05 Ha/Bohr²。它是解释采样和接口的教学模型。真实分子、材料或界面应改用合适的电子结构方法或力场，并重新检查采样条件。配套运行与比较脚本对应本文固定对照；正式项目应按体系、温度与目标量调整输入和分析。

两个体系目标温度分别为300 K和600 K，单珠子、NVT、Langevin时间常数10 fs，步长0.5 fs、4000步，共2 ps。这里展示经典单珠子温度交换；路径积分珠数教程另行说明量子环聚合物。

3 i-PI输入：定义两个温度与交换#

下面是操作包的实际控制模板inputs/input.xml。启动脚本会为当前工作目录生成独立的本地连接名称：

<simulation verbosity="low">
 <output prefix="remd">
  <properties stride="1" filename="out">[step,time{femtosecond},potential,kinetic_md,temperature{kelvin},ensemble_temperature{kelvin},conserved]</properties>
  <trajectory filename="pos" stride="1">positions{angstrom}</trajectory>
  <trajectory filename="force" stride="1">forces</trajectory>
  <checkpoint stride="100"/>
 </output>
 <total_steps>4000</total_steps>
 <prng><seed>20261109</seed></prng>
 <ffsocket name="cp2k" mode="unix" matching="any" pbc="false">
  <address>ys438-20f20aa640f28b80</address><latency>0.00001</latency><timeout>60</timeout>
 </ffsocket>
 <system_template>
  <labels>[IDX,TEMP]</labels>
  <instance>[0,300]</instance><instance>[1,600]</instance>
  <template>
   <system prefix="replica-IDX">
    <initialize nbeads="1">
     <file mode="xyz" units="angstrom">init.xyz</file>
     <masses mode="manual" units="dalton">[1.0]</masses>
     <velocities mode="thermal" units="kelvin">TEMP</velocities>
    </initialize>
    <forces><force forcefield="cp2k"/></forces>
    <motion mode="dynamics">
     <fixcom>false</fixcom>
     <dynamics mode="nvt">
      <timestep units="femtosecond">0.5</timestep>
      <thermostat mode="langevin"><tau units="femtosecond">10</tau></thermostat>
     </dynamics>
    </motion>
    <ensemble><temperature units="kelvin">TEMP</temperature></ensemble>
   </system>
  </template>
 </system_template>
 <smotion mode="remd"><remd><stride>20</stride><krescale>true</krescale><swapfile>exchange.dat</swapfile></remd></smotion>
</simulation>

system_template定义两个体系；TEMP分别替换为300、600。fixcom=false保留外部势中粒子的运动自由度，初始位置文件显式使用Å。单粒子初态为：

1
# CELL(abcABC): 20 20 20 90 90 90
H 0.1 0.1 0.1

remd/stride=20表示每个MD步按1/20概率尝试交换，本例只有一对温度；它不表示严格每20步交换一次。krescale=true在交换温度后同步缩放动量。

4 CP2K输入：作为力计算客户端#

两个客户端分别在results/client-0/和results/client-1/运行。它们采用相同势能方法，为i-PI请求的坐标计算能量和力。以下为第一个模板的完整输入：

&GLOBAL
 PROJECT cp2k-client-0
 RUN_TYPE DRIVER
 PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
 &DRIVER
  HOST ys438-20f20aa640f28b80
  UNIX TRUE
  SLEEP_TIME 0.00001
 &END DRIVER
&END MOTION
&FORCE_EVAL
 METHOD FIST
 &MM
  &FORCEFIELD
   DO_ELECTROSTATICS FALSE
   &CHARGE
    ATOM H
    CHARGE 0
   &END CHARGE
   &NONBONDED
    &LENNARD-JONES
     ATOMS H H
     EPSILON 0.0
     SIGMA [angstrom] 1
     RCUT [angstrom] 4
    &END LENNARD-JONES
   &END NONBONDED
  &END FORCEFIELD
  &POISSON
   PERIODIC NONE
   &EWALD
    EWALD_TYPE NONE
   &END EWALD
  &END POISSON
 &END MM
 &EXTERNAL_POTENTIAL
  ATOMS_LIST 1
  FUNCTION 0.5*K*(X^2+Y^2+Z^2)
  PARAMETERS K
  VALUES 0.05
 &END EXTERNAL_POTENTIAL
 &SUBSYS
  &CELL
   ABC 20 20 20
   PERIODIC NONE
  &END CELL
  &TOPOLOGY
   CONNECTIVITY OFF
  &END TOPOLOGY
  &KIND H
   ELEMENT H
   MASS 1.0
  &END KIND
  &COORD
   H 0.1 0.1 0.1
  &END COORD
 &END SUBSYS
&END FORCE_EVAL

关键是RUN_TYPE DRIVER、MOTION/DRIVER的连接名和UNIX TRUE。本次CP2K2026.2客户端的UNIX命名空间为/tmp/qiskit_，配套i-PI启动时使用-S /tmp/qiskit_匹配；操作包已设置好这一选项。自己组织输入时，客户端名称与服务端地址必须一致。这里使用本机UNIX连接。自己的周期体系应在i-PI初态中明确晶胞单位，CP2K客户端使用服务端传来的晶胞和坐标。

5 实际发生了哪些交换#

source env.sh
"$YS_PY" - <<'PYCODE'
import json
from pathlib import Path
v = json.loads(Path('results/VALIDATION.json').read_text())
print('实际尝试:', v['attempts'])
print('接受:', v['accepted'], '拒绝:', v['rejected'])
print('交换记录: results/remd.exchange.dat')
PYCODE

本文参考轨迹中实际尝试185次，接受99次。两种温度、相同势能、缩放动量时，交换概率可写为min(1, exp[(βᵢ−βⱼ)(Uᵢ−Uⱼ)])。这是基于真实当前势能的判定，交换不必每次接受。

i-PI3.3.1默认对多个体系并行推进，随机数的调用顺序会随线程执行顺序变化。因此，即使保留种子，重新运行的坐标轨迹和接受次数也可能与参考图不同。应检查本次真实能量、力、温度访问和交换记录；若希望按固定调用顺序操作，可按i-PI文档研究threading=false的设置。

6 下一步：按温度分析#

replica-0_remd.out代表物理walker 0，里面的目标温度会发生改变。固定温度的平均能量或结构分布应先按真实交换历史重新整理，见按温度整理副本交换结果。2 ps用于展示操作，正式研究还需检查充分混合、往返次数、自相关与统计稳定性。

参考：i-PI温度交换、i-PI输入字段、CP2K DRIVER、本版本UNIX连接实现。