PLUMED高斯核密度：带宽、归一化和积分

核密度估计在每个样本位置放一个平滑核，再叠加归一化。本篇用原生KDE直接处理与离散直方图相同的32个距离值，便于比较两种表示。

1. 下载和实际密度#

下载本篇完整输入、原生网格和绘图资料包，保存为example.zip：

python3 -m zipfile -e example.zip .
cd plumed-histogram-grid-examples
bash results.sh
bash read.sh
bash plot.sh

这些命令先读取已完成的实际结果；原生输出在reference-results/，CSV在data/，图在figures/。默认使用集群ASE环境Python，可用YS_ANALYSIS_PYTHON指定含NumPy、ASE和Matplotlib的解释器。

这是32个人工数值的工具例，不称为真实MD平衡概率密度。

2. 完整输入#

cat inputs/gaussian-kde-explicit-height/plumed.dat

UNITS LENGTH=nm ENERGY=kj/mol TIME=ps
values: CONSTANT VALUES=0.05000000000000,0.06290322580645,0.07580645161290,0.08870967741935,0.10161290322581,0.11451612903226,0.12741935483871,0.14032258064516,0.15322580645161,0.16612903225806,0.17903225806452,0.19193548387097,0.20483870967742,0.21774193548387,0.23064516129032,0.24354838709677,0.25645161290323,0.26935483870968,0.28225806451613,0.29516129032258,0.30806451612903,0.32096774193548,0.33387096774194,0.34677419354839,0.35967741935484,0.37258064516129,0.38548387096774,0.39838709677419,0.41129032258065,0.42419354838710,0.43709677419355,0.45000000000000
amplitude: CONSTANT VALUE=26.59615202676218
kde: KDE ARG=values GRID_MIN=0 GRID_MAX=0.5 GRID_BIN=100 KERNEL=GAUSSIAN BANDWIDTH=0.015 CUTOFF=1000 HEIGHTS=amplitude
density: CUSTOM ARG=kde FUNC=x/32 PERIODIC=NO
area: INTEGRATE_GRID ARG=density PERIODIC=NO
PRINT ARG=area FILE=COLVAR STRIDE=1 UPDATE_FROM=1 FMT=%.14g
DUMPGRID GRID=density FILE=density.dat STRIDE=0 FMT=%.14g

values是一个32维常量向量。这里直接对完整样本快照做KDE，反复打印不会向样本集合额外添加数据；它不同于HISTOGRAM沿轨迹逐帧累计。

3. 带宽与核高度#

本例σ=0.015 nm，单个核为exp(−(x−xᵢ)²/(2σ²))/(sqrt(2π)σ)；本例32个核求和后再除32。

HEIGHTS=amplitude显式给足有效数字的核高度，使本例与独立公式的对照不受快捷入口中归一化常数舍入影响。HEIGHTS是峰高，不能直接把32个样本的峰高都设为1后仍称密度已经归一化。

更换带宽时同步重算高度，例如0.02 nm：

python3 - <<'PYCODE'
import math
sigma = 0.02
print(1 / (math.sqrt(2 * math.pi) * sigma))
PYCODE

同时修改BANDWIDTH及相应样本数量；此包的read.py是固定教学条件的对照脚本，更换条件后也要同步分析参数。带宽改变平滑尺度，小带宽并不自动使统计结果更准确。

CUTOFF=1000在这个小例中使全部有限网格点均进入高斯计算，便于解析对照；大型多维网格应考虑实际计算量和适当截断。

4. 运行与提取#

在新解压目录执行短串行driver：

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

环境使用CP2K2026.2依赖中的PLUMED2.10.1，run.py已设置相应程序和动态库路径。这些是坐标/网格处理，不做电子结构或MD积分。

head reference-results/gaussian-kde-explicit-height/density.dat
cat reference-results/gaussian-kde-explicit-height/COLVAR
head data/PLUMED-Gaussian-KDE-grid.csv

区间0到0.5 nm、100个间隔，共101个节点。密度单位nm⁻¹，导数单位nm⁻²。本例原生密度与独立公式最大差约5.11×10⁻¹⁴，原生导数与解析导数最大差约4.83×10⁻¹³。

5. 归一化与数值积分#

本例网格积分为0.9999859429；同一区间的解析高斯积分为0.9999723119，相差约1.36×10⁻⁵。高斯尾部超出有限区间，网格积分又有离散误差，所以不能要求打印结果数学上精确等于1。

实际轨迹还要考虑样本相关性、平衡段和权重，不能仅因KDE平滑就称结果收敛。查网格文件和插值的实际操作，见参考网格与查询。