LAMMPS Drude位点与Thole屏蔽

Drude模型用带电辅助粒子相对核心的位移表示极化。两个位点靠近时，短程相互作用的处理很重要。本篇从两个固定几何位点出发，演示如何叠加Thole屏蔽并读取能量。

下载五篇共用的完整输入、科学数据与绘图包。先按照多电子波包篇的共用运行步骤提交一次，后续四篇读取同一批结果，不需要重复提交。

1. 准备两个中性位点#

cat inputs/thole-0p8/system.data
cat inputs/bare-0p8/input.lammps

完整结构：

Owned core-Drude model

4 atoms
2 bonds
2 atom types
1 bond types

-5 15 xlo xhi
-5 15 ylo yhi
-5 15 zlo zhi

Masses

1 11.6
2 0.4

Atoms # full

1 1 1 0.5 0 0 0
2 1 2 -0.5 0.050000000000000003 0 0
3 2 1 0.5 0.80000000000000004 0 0
4 2 2 -0.5 0.85000000000000009 0 0

Bonds

1 1 1 2
2 1 3 4

每个位点核心+0.5e、Drude−0.5e，质量分别11.6u/0.4u。各Drude沿+x偏移0.05Å，谐键K=500kcal·mol⁻¹·Å⁻²，r0=0。LAMMPS谐键用 K(r−r0)²，因此两个键合计2.5kcal/mol；不要再加1/2。

这两个位移是输入给定，尚未自洽最小化，不把它称为电场下已求出的诱导响应。单个位点在电场下的最小化可参考Drude谐响应。

2. 完整Thole输入#

units real
atom_style full
atom_modify map array
boundary p p p
read_data system.data
bond_style harmonic
bond_coeff 1 500 0
special_bonds lj/coul 0 0 0
fix roles all drude C D
pair_style hybrid/overlay coul/cut 8.0 thole 2.6 8.0
pair_coeff * * coul/cut
pair_coeff * * thole 0.0830159275
comm_modify vel yes
thermo 1
thermo_style custom step pe ebond ecoul evdwl
thermo_modify norm no format float %.17g
run 0
write_dump all custom final.dump id mol type q x y z fx fy fz modify sort id format float %.17g

fix roles all drude C D依次指定核心类型、Drude类型。thole与Coulomb通过 hybrid/overlay叠加；它不是独立替代全部电荷相互作用的完整分子力场。

本例polarizability按当前real单位换算为 332.06371×0.5²/(2×500)=0.0830159275ų，阻尼参数2.6。这里只演示固定参数的模型行为，不推荐把这套质量、电荷、键强或极化率直接当作某种真实原子的参数。

3. 三个距离的实际数据#

位点间能是总势能减去两个内部谐键能。裸电荷条件还由四个跨位点 332.06371 qi qj/rij 项独立求和。图中右侧为0.8Å条件的实际坐标。

4. 不要漏掉键力#

cat inputs/thole-0p8/final.dump
cat inputs/thole-FD-minus/input.lammps
cat inputs/thole-FD-plus/input.lammps

对0.8Å的裸电荷和Thole两组，都实际将第一个核心x改变±1e−5Å，用总能中心差分对照fx。核心移动也改变本位点的谐键长度，因此必须使用含键能的总势能，不能只对上表的非键能差分。

已有原生结果时，只重新提取和画图：

module load sisl/0.16.4
python analyze.py
python plot.py

这一步读取实际结果，不重新运行原生模拟。

参考LAMMPS Thole用法。如要做有限温度可极化轨迹，继续双温控与两种温度。