GROMACS 外电场与偶极响应

外加电场可以改变带电原子的受力，并通过分子部分电荷影响取向。本篇用同一SPC/E水盒比较零场、+z与−z静态电场，先弄清输入、单位、偶极和原子力，再讨论能够从短轨迹中读到什么。

下载完整输入、真实轨迹与绘图资料包，约147 MiB。资料包含32个短算例及6次同构型电场力对照。解压进入资料包目录，先查看已有结果：

bash view.sh

四张轨迹数据图位于 analysis/figures，电场力图位于 field-analysis/figures；每个条件的CSV和数值汇总位于 analysis/条件名。脚本调用集群已有MDAnalysis2.9.0、NumPy、Matplotlib环境，画已有数据时不提交计算。

本例使用348个SPC/E水分子，共1044个原子，初始坐标来自已有水盒的NVT平衡末帧。新计算重新生成速度，300 K、固定盒子NVT、步长1 fs、长度20 ps，每0.1 ps保存一帧，共201帧完整位置/速度TRR。两种速度种子13579、24680分别在P100和V100 16G计算。不同平台重复相同种子，不当作四个独立初始种子。

GROMACS2026.4 thread-MPI，GNU13.3、CUDA12.8、FFTW3.3.10。短程非键和PME使用GPU，键相互作用及坐标更新使用CPU；Colvars集体变量的计算不因此变成GPU内核。

图中箭头依据实际分子部分电荷和坐标计算，再统一缩放显示方向；箭头长度不是分子尺寸。

1 准备与提交#

python3 prepare.py
cd cases/p100-field-plus-s13579
cat md.mdp
sbatch job.slurm
squeue -u "$USER"

零场目录为 p100-zero-s13579，负场目录为 p100-field-minus-s13579；V100与第二种子也已配套。流程：同一平衡坐标 → 重新生成初速度 → 指定场方向 → 20 ps NVT → 原生偶极提取 → 与零场、反向场对照。

模板每项申请1张对应GPU和1个CPU核心，使用核时分区；套餐用户按账户规则修改分区。不要在登录节点直接运行完整MD。已有 cases 或结果目录时，重新解压到新目录准备，避免覆盖。

任务实际执行的主命令如下，供已分配资源中的手动操作参考：

gmx grompp -f md.mdp -c initial.gro -p topol.top -o md.tpr   -po md-processed.mdp -pp md-processed.top
gmx mdrun -deffnm md -ntmpi 1 -ntomp 1 -nb gpu -pme gpu   -bonded cpu -update cpu

拓扑完整分子表如下；原子电荷、刚性水几何由所引用的SPC/E模型给出：

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/spce.itp"

[ system ]
SPC/E liquid transport tutorial in water

[ molecules ]
SOL               348

2 电场参数和完整MDP#

integrator = md
dt = .001
nsteps = 20000
cutoff-scheme = Verlet
nstlist = 20
coulombtype = PME
rcoulomb = .9
rvdw = .9
DispCorr = EnerPres
fourierspacing = .12
pme-order = 4
pbc = xyz
constraints = none
nstenergy = 100
nstcalcenergy = 100
nstlog = 1000
nstxout = 100
nstvout = 100
nstfout = 0
nstxout-compressed = 0
tcoupl = v-rescale
tc-grps = System
tau-t = .5
ref-t = 300
pcoupl = no
comm-mode = Linear
nstcomm = 100
continuation = no
gen-vel = yes
gen-temp = 300
gen-seed = 13579
ld-seed = 13680
electric-field-z = .05 0 0 0

electric-field-z依次为 E0 omega t0 sigma，单位分别是V/nm、ps⁻¹、ps、ps。本例 .05 0 0 0 是+z方向0.05 V/nm静态场；改成 -.05 0 0 0 为反向场。零场条件不添加这行。

omega = 0表示不振荡，sigma = 0表示不加高斯脉冲包络。不要把V/nm当作V/Å，也不要把角频率直接当成Hz。

3 带电原子受到什么力#

F_z [kJ/(mol·nm)] = q [e] * E_z [V/nm] * 96.48533212

SPC/E中氧电荷−0.8476e，每颗氢+0.4238e。+z电场下，氧的额外力指向−z，氢指向+z；电中性的单个水分子净附加力为0，但仍可受到取向力矩。

力单位为kJ/(mol·nm)。附件还在完全相同坐标上，用两种构建分别做零场、正场、负场单帧CPU重算，直接比较“有场总力−零场总力”。这是CPU力对照，不是GPU速度比较。

全部1044个原子的最大分量差小于0.000182 kJ/(mol·nm)。图展示前18个原子以便看清氧、氢交替的力方向，完整原始力随附件提供。

本版本此外场通过附加力施加；同构型对照中，原生 Potential 数值保持相同。不能把有场任务的原生Potential直接当成包含完整外场势能的热力学自由能。场下轨迹的能量也会受坐标变化、外场做功及恒温器影响。

4 任务结束后提取结果#

在当前算例目录执行：

module load gromacs/2026.4-p100-threadmpi
python3 ../../export_native.py
ls native-analysis

V100更换为对应模块。这个步骤读取已有轨迹和能量，输出到新建的 native-analysis，不会覆盖同名分析目录。全部参考结果及绘图脚本可在资料包根目录运行 bash view.sh 查看。

原生偶极提取可单独执行：

printf 'SOL
' | gmx dipoles -s md.tpr -f md.trr   -o dipole.xvg -eps dielectric.xvg -temp 300

附件的 native-analysis/dipole.xvg 列为时间、Mx、My、Mz、总偶极模，偶极单位Debye。程序对分子部分电荷求和；处理自己的体系时，应先保证分子完整和正确的电荷定义。带自由离子的体系还涉及跨周期边界与导电贡献，不能照搬中性水的处理。

在当前条件目录画自己的新结果：

module load ase/3.29.0
python - <<'PY'
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
d=np.loadtxt("native-analysis/dipole.xvg",comments=["#","@"],ndmin=2)
plt.plot(d[:,0],d[:,3])
plt.xlabel("Time (ps)");plt.ylabel("Total water dipole Mz (Debye)")
plt.tight_layout();plt.savefig("new-field-dipole.png",dpi=180)
PY

这段程序直接读取刚生成的原生数据；PNG保存在当前目录。XVG以 # 和 @ 开头的行是说明/列名，不能当作数值读入。

可选：重新做同构型力对照#

在已分配的CPU资源中使用1个线程。下面读取附件初始构型，分别执行零场、正场、负场单帧原生重算；不延长MD轨迹。先回到资料包根目录：

module load gromacs/2026.4-p100-threadmpi
PY=/public/toolkit/miniconda3/envs/qvasp-2.26/bin/python
mkdir field-controls-new
"$PY" run_field_controls.py --platform p100 --output field-controls-new/p100
python3 export_field_controls.py --controls field-controls-new/p100
"$PY" plot_field.py --controls field-controls-new --platform p100 --output field-check-new

field-check-new/figures得到原生力和qE对照图，SUMMARY.json给出全部原子的最大分量差。已有同名目录时使用新目录，不覆盖参考结果。V100对应更换模块和platform参数。

5 零场与正负场实际结果#

本篇读左列；同图中间和右列是非平衡余弦加速的配套数据。

正场结果的Mz均值为正，负场为负；零场短窗口仍有明显涨落。这里看的是348个水分子的总偶极，不是每分子平均，也不是介电常数。20 ps足以演示输入与输出，但不足以给出可靠介电响应、线性响应区间或误差条。

要做材料结论，需检查场强依赖、平衡时间、采样长度、体系尺寸、恒温方式和不同初始种子，并针对目标定义极化或介电响应。本页0.05 V/nm是操作示例参数，不是所有体系的推荐场强。

完整参数见GROMACS电场MDP选项，偶极分析见原生dipoles说明。