PLUMED平滑配位数与接触图

配位数作为集体变量时，通常需要连续可导的“平滑计数”，而不是距离一过阈值就突然从1变成0。本篇用32帧人工构造的四粒子几何演示输入与数据读取，没有进行物理MD或自由能采样。

1. 下载并看结构#

下载完整输入、原生输出、CSV与绘图脚本，保存为example.zip：

python3 -m zipfile -e example.zip .
cd plumed-contacts-coordination
module load ase/3.29.0
bash read.sh
bash plot.sh

这一步读取本包保存的数据，在analysis-read/生成结果表、在figures/生成图片。新运行后处理时使用bash run.sh，输出写入新建的work-new/；它不运行VASP，也不提交新作业。

这里调用CP2K2026.2依赖目录中的PLUMED2.10.1，不是把公共独立PLUMED2.9.x模块当成同一个程序。run.sh只运行PLUMED driver，已设置对应动态库路径。

长度为nm，能量为kJ/mol。粒子编号1为中心，分别与2、3、4计算接触；这些编号是1起始的原子序号。

2. 核心输入#

显式使用同一个切换函数：

UNITS LENGTH=nm ENERGY=kj/mol TIME=ps
cn: COORDINATION GROUPA=1 GROUPB=2-4 SWITCH={RATIONAL R_0=0.3 NN=6 MM=12} NOPBC
cp: COORDINATION GROUPA=1,2 GROUPB=3,4 SWITCH={RATIONAL R_0=0.3 NN=6 MM=12} PAIR NOPBC
cm: CONTACTMAP ATOMS1=1,2 ATOMS2=1,3 ATOMS3=1,4 SWITCH={RATIONAL R_0=0.3 NN=6 MM=12} NOPBC
PRINT ARG=cn,cp,cm.* FILE=COLVAR STRIDE=1

这组参数的接触函数可写成1/(1+(r/0.3)^6)，距离0.3 nm时为0.5，不是“0.3 nm以内全为1、以外全为0”。cn把1–2、1–3、1–4三个接触相加；PAIR只把两组按顺序配对，所以cp计算1–3与2–4。

3. 实际运行和读取#

在完整资料包中：

bash run.sh
module load ase/3.29.0
python -B read.py work-new analysis-new
python -B plot.py analysis-new figures-new
head work-new/explicit-switch/COLVAR

work-new/保留简写参数入口，work-new/explicit-switch/保留显式SWITCH入口。以首帧为例，三个接触约为0.83986965、0.41221160、0.17939485；显式SWITCH的配位数约为1.43147610，与三者之和一致。

4. 为什么两种写法末位不一样#

本机2.10.1中，直接写R_0=0.3 NN=6 MM=12的简写入口会自动引入有限截断并做拉伸；本例显式SWITCH字符串未给D_MAX，没有同样的有限截断拉伸。

32帧中两种计数的最大差约1.62×10⁻⁵。这来自默认函数定义，不是接触图漏了原子。需要把配位数与接触分量逐项对应时，像本篇一样给两者写同一个完整SWITCH，并检查实际输出；不能仅因R_0相同就认为所有默认设置相同。

5. 邻居表的使用范围#

完整输入另有：

nl: COORDINATION GROUPA=1 GROUPB=2-4 SWITCH={RATIONAL R_0=0.3 NN=6 MM=12} NOPBC NLIST NL_CUTOFF=1.0 NL_STRIDE=4

本例所有相关距离均小于1 nm，因此与完整求和一致。实际大体系需要根据切换函数、移动速度和更新间隔设置邻居表，不能把本例1 nm当成通用截断。

离线driver还传入--box 2,2,2满足此版本邻居表的盒元数据要求；本例CV均明确NOPBC，不把这条盒声明当成已开展周期MD。

PLUMED配位数说明给出分组与切换函数定义。参考接触图和偏置力见下一篇。