PLUMED：Karplus关系与J耦合

本篇演示 JCOUPLING TYPE=CUSTOM：从四原子二面角计算经验 Karplus 关系，并读取 J 耦合和约束力。经验参数依赖具体核对和化学环境，本例参数仅用于清楚展示操作。

1. 下载并提取#

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

python3 -m zipfile -e example.zip .
cd plumed-FRET-RDC-J
module load ase/3.29.0
bash read.sh
bash plot.sh

本包用 PLUMED 2.10.1，调用集群 CP2K 2026.2 的独立依赖程序，不运行 CP2K 电子结构计算。输入是 32 个四粒子几何帧，用来演示可观测量和偏置力，不是蛋白质 MD 或实验数据。

默认读取 reference-results/。若需要重算，在根目录运行 bash submit.sh：CPU 模板申请 56 核整节点，PLUMED 使用 1 线程。作业结束后执行 bash read.sh work-new analysis-read，再 bash plot.sh。一次作业同时提取 FRET、两组 RDC 和一组 J 耦合。

2. 写出参数和原子顺序#

UNITS LENGTH=nm ENERGY=kj/mol TIME=ps
j: JCOUPLING TYPE=CUSTOM A=6.4 B=-1.4 C=1.9 SHIFT=0.3 ATOMS1=1,2,3,4 NOPBC
t: TORSION ATOMS=1,2,3,4 NOPBC
PRINT ARG=j.*,t FILE=Jcoupling.dat STRIDE=1

关系是 J=A·cos²(θ+SHIFT)+B·cos(θ+SHIFT)+C。这里 SHIFT=0.3 使用 弧度，不是 0.3 度；原生 TORSION 输出也用弧度，画图脚本才转换为度。

四原子的顺序必须与参数定义一致。把 ψ、φ 或其他二面角混用，即使程序正常结束，也不代表算对了目标核对。可选的内置类型及自定义关系参见 官方 JCOUPLING 文档。

3. 运行和定位输出#

cat inputs/plumed.dat
bash submit.sh

计算后：

bash read.sh work-new analysis-read
bash plot.sh
head -6 reference-results/COLVAR

本包只有一组四原子定义，输出是 j.j-0。实际 32 帧的 J 值约 2.684595–3.167285 Hz。图中的整条灰线是所选经验公式的解析曲线，蓝点才是本次实际几何的原生输出；不要把整条曲线误认为跑过了所有二面角。

4. 把J耦合接入偏置#

本包实际输入还包含：

bj: RESTRAINT ARG=j.j-0 AT=4 KAPPA=0.3

目标 4 Hz，KAPPA 对应 kJ/(mol·Hz²)。driver 输出 FRET、RDC 和 J 三个约束产生的总偏置力，本次逐一核对了 32帧×4粒子×3方向=384 个分量，与独立几何公式的中心差分最大差约 1.19×10⁻⁹ kJ/(mol·nm)。

head -8 analysis-read/PLUMED-spectroscopy-bias-forces.csv

这个比较说明本例输入下的量和力传递一致。它不能替代经验参数对真实分子和实验的适用性判断。

5. 换成自己的体系#

先确定核对、对应二面角和参数来源，再核对原子编号、角度正负与周期边界；用已知构型检查几条输出，最后按需要做时间平均或带权平均。经验 J 耦合与距离效率一样具有非线性，不能直接用平均结构代替对轨迹逐帧计算。

若只需计算可观测量，不要加入本篇的 RESTRAINT 行；若需要实验约束采样，再设置有依据的目标、强度和统计方式。