LAMMPS Drude振子与线性极化

Drude模型把一个可极化位点拆成核与辅助带电粒子，通过位移产生诱导偶极。本例先用可精确计算的单个谐振子理解它的结构、单位和电场响应，再考虑复杂分子力场。

下载三篇共用的完整输入、实算数组、原生科学输出和绘图包。按共用运行步骤运行一次。这里没有Thole屏蔽或双温热化，不是某种溶剂的现成力场参数。

1. 完整两粒子结构#

核质量11.6u、电荷+0.5e，Drude粒子质量0.4u、电荷−0.5e，总质量12u、总电荷0。类型1/2分别是核与Drude，不代表两种化学元素。

cat inputs/drude-one/oscillator.data
cat inputs/drude-one/input.lammps

Harmonic neutral core-Drude model

2 atoms
1 bonds
2 atom types
1 bond types

-5 5 xlo xhi
-5 5 ylo yhi
-5 5 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.001 0 0

Bonds

1 1 1 2

2. 完整最小化输入#

units real
atom_style full
atom_modify map array
boundary f f f
read_data oscillator.data
pair_style zero 6.0
pair_coeff * *
bond_style harmonic
bond_coeff 1 500.0 0.0
special_bonds lj/coul 0 0 0
fix polar all drude C D
fix field all efield 1 0 0
fix_modify field energy yes
compute dipole all dipole
thermo 1
thermo_style custom step pe c_dipole[1] c_dipole[2] c_dipole[3]
thermo_modify norm no format float %.17g
min_style cg
minimize 0 1e-12 200 1000
write_dump all custom final.dump id mol type q x y z fx fy fz modify sort id format float %.17g
write_data final.data nocoeff

fix polar all drude C D标记两类粒子。bond_style harmonic把它们连起来，K=500 kcal/(mol·Å²)、平衡长度0；LAMMPS的谐键势在这里是 K×r²，已经把常见的1/2吸收到K里，因此求位移时分母是 2K。

fix field ... efield 1 0 0施加x向1 V/Å的均匀电场。fix_modify field energy yes把电场势能计入最小化；否则只看力与只看势能的终止条件可能混淆。非键相互作用为零，核与辅助粒子之间由谐键代表该模型内部的极化约束。

3. 位移、偶极和单位#

定义相对位移 d=xDrude−xCore，解析结果为 d=qDrude×23.060549×E/(2K)。23.060549是本版本real单位中电场到力的换算系数；电荷用e，电场V/Å，K用kcal/(mol·Å²)，得到d为Å。

中性位点偶极 p=Σqᵢxᵢ=qDrude×d，单位e·Å。负Drude粒子向电场反方向移动，但诱导偶极沿电场方向；移动整个中性位点不会改变偶极。

四个条件都实际进行了最小化。正负电场下响应反向，零场位移趋于0，偶极与电场成正比。脚本用完整电荷和坐标重算偶极，也比较 K|d|²−23.060549×E×p 与原生总势能。

cat inputs/drude-one/final.dump
module load sisl/0.16.4
python analyze.py
python plot.py

复杂分子模型还需要可靠的电荷、极化率、阻尼屏蔽、短程势、核心/辅助粒子质量划分和合适的控温。这里的静态线性关系不证明有限温度Drude动力学、分子间屏蔽或实际介电常数；不要把这套演示参数直接用于材料预测。

参考：Drude粒子标记、可极化模型说明。