• 科研
• 科学计算
• 高性能计算

pair_style spin/exchange与pair_style spin/exchange/biquadratic是 LAMMPS SPIN 包中用于描述原子磁自旋之间交换相互作用（exchange interaction）的核心对势风格。本文以官方文档 pair_spin_exchange.rst 为主体，系统讲解两种风格的语法、Heisenberg 型哈密顿量与 Bethe-Slater 型径向参数化函数、磁力矩与机械力的推导、pair_coeff参数顺序、offset能量偏移选项，并结合仓库源码与examples/SPIN中的铁、镍、钴实际输入文件，给出从参数拟合到可运行输入脚本的完整实战方案。读完本文，你将掌握如何在atom_style spin体系中正确配置磁交换相互作用、理解其底层实现原理，并能为自旋-晶格（spin-lattice）动力学与磁性能计算（如状态方程）选择正确的参数。

一、命令语法与快速上手示例

两种风格的基本语法完全一致，仅有一个全局截断距离参数（单位采用 metal 单位制的距离单位，即 Å）：

pair_style spin/exchange cutoff pair_style spin/exchange/biquadratic cutoff

其中cutoff为全局截断距离（metal 单位制下的距离单位）。

官方文档给出的典型示例：

pair_style spin/exchange 4.0 pair_coeff * * exchange 4.0 0.0446928 0.003496 1.4885 pair_coeff 1 2 exchange 6.0 -0.01575 0.0 1.965 offset yes pair_style spin/exchange/biquadratic 4.0 pair_coeff * * biquadratic 4.0 0.05 0.03 1.48 0.05 0.03 1.48 offset no pair_coeff 1 2 biquadratic 6.0 -0.01 0.0 1.9 0.0 0.1 19

注意pair_coeff中* *之后的第一个字符串是子风格名：spin/exchange使用exchange，spin/exchange/biquadratic使用biquadratic。从源码src/SPIN/pair_spin_exchange.cpp的coeff()函数可以看到，该子风格名是强制校验的：if (strcmp(arg[2],"exchange") != 0) error->all(FLERR,"Incorrect args for pair coefficients" ...)，且参数个数必须为 7 或 9（9 个参数对应附加offset关键字）源码见。biquadratic风格则校验arg[2]必须为"biquadratic"，参数个数为 10 或 12 源码见。

二、物理模型：两种交换相互作用哈密顿量

2.1 经典 Heisenberg 交换项（spin/exchange）

spin/exchange风格计算一对磁自旋之间的交换相互作用，哈密顿量为：

$$ H_{ex} = -\sum_{i,j}^N J_{ij} (r_{ij}) ,\vec{s}_i \cdot \vec{s}_j $$

其中 $\vec{s}_i$ 与 $\vec{s}j$ 是两个粒子（通常是原子）磁自旋的单位向量，$r{ij} = |\vec{r}_i - \vec{r}j|$ 是两粒子间的原子间距，求和覆盖最近邻对。$J(r{ij})$ 是一个描述不同近邻壳层上交换相互作用强度与符号的函数——$J>0$ 对应铁磁性（自旋平行能量更低），$J<0$ 对应反铁磁性倾向。

2.2 含双二次项的扩展（spin/exchange/biquadratic）

spin/exchange/biquadratic风格在 Heisenberg 项之外额外引入双二次（biquadratic）交换项：

$$ H_{bi} = -\sum_{i, j}^{N} {J}{ij} \left(r{ij} \right), \vec{s}{i}\cdot \vec{s}{j} -\sum_{i, j}^{N} {K}{ij} \left(r{ij} \right), \left(\vec{s}{i}\cdot \vec{s}{j}\right)^2 $$

其中 $\vec{s}i$、$\vec{s}j$、$r{ij}$ 与 $J(r{ij})$ 的定义与上文一致，$K(r_{ij})$ 是第二个函数，定义双二次项的强度与符号。双二次项对自旋取向的非线性响应能够描述更复杂的磁序（如某些多铁性材料与自旋玻璃行为），examples/SPIN目录中的 BFO（铁酸铋）算例即使用多种自旋相互作用叠加建模多铁性材料。

三、Bethe-Slater 型径向参数化函数

两种相互作用中 $J(r_{ij})$ 与 $K(r_{ij})$ 随原子间距的变化都由同一个函数定义：

$$ f\left( r_{ij} \right) = 4 a \left( \frac{r_{ij}}{d} \right)^2 \left( 1 - b \left( \frac{r_{ij}}{d} \right)^2 \right) e^{-\left( \frac{r_{ij}}{d} \right)^2 }\Theta (R_c - r_{ij}) $$

其中 $a$、$b$、$d$ 是三个在相应pair_coeff命令中定义的常系数，$R_c$ 是与该对相互作用关联的截断半径，$\Theta(R_c - r_{ij})$ 为阶跃函数，保证 $r_{ij} \ge R_c$ 时相互作用严格为零。该函数本质上是著名的 Bethe-Slater 曲线的解析拟合形式：因子 $e^{-(r_{ij}/d)^2}$ 提供平滑衰减，$1 - b(r_{ij}/d)^2$ 提供符号翻转能力，使同一函数既能描述铁磁近邻壳层也能描述反铁磁近邻壳层。

系数 $a$、$b$、$d$ 需要通过拟合确定，使上述函数与所考虑的 $N$ 个近邻壳层上的交换相互作用数值相匹配。对于spin/exchange/biquadratic风格，需要拟合六个系数（$J(r_{ij})$ 三个、$K(r_{ij})$ 三个）。

3.1 从源码看函数的具体求值

src/SPIN/pair_spin_exchange.cpp中的compute_exchange()函数给出了该函数在代码中的精确实现 源码见：

ra = rsq/J3[itype][jtype]/J3[itype][jtype]; Jex = 4.0*J1_mag[itype][jtype]*ra; Jex *= (1.0-J2[itype][jtype]*ra); Jex *= exp(-ra);

这里rsq/J3/J3对应 $(r_{ij}/d)^2$，J1_mag即系数 $a$（换算为磁力矩时除以约化普朗克常数 $\hbar$，对应构造函数中的J1_mag[i][j] = j1/hbar），J2即系数 $b$，J3即系数 $d$。三行代码与文档公式逐项对应：4·a·(r/d)²·(1 − b·(r/d)²)·exp(−(r/d)²)。biquadratic风格的compute_exchange()则对 $J$ 与 $K$ 各求一次该函数，并将自旋点积及其平方项组合进磁力矩 源码见。

3.2 使用仓库脚本进行参数拟合

examples/SPIN目录为铁（bcc）、镍（fcc）、钴（fcc/hcp）等材料提供了完整的交换参数拟合脚本。以镍为例，exchange_fit.py 使用 SciPy 的curve_fit对 Bethe-Slater 函数：

def func(x,a,b,c): return 4*a*((x/c)**2)*(1-b*(x/c)**2)*np.exp(-(x/c)**2)

进行拟合，数据文件 exchange_fcc_ni.dat 提供各近邻壳层距离与对应的 $J$ 值（如2.492 0.0028027、3.524 0.0000816等），拟合结果即为pair_coeff中的 $(a,b,d)$。examples/SPIN/README明确指出这些算例的目的在于演示使用方法，真实材料模拟需要更精细的参数调优。

四、磁力矩与机械力：自旋-晶格耦合的桥梁

由交换相互作用，每个自旋 $i$ 会受到磁力矩 $\vec{\omega}{i}$，其关联原子在自旋-晶格计算（见 fix nve/spin）中还会受到力 $\vec{F}{i}$：

$$ \vec{\omega}{i} = \frac{1}{\hbar} \sum{j}^{Neighb} {J} \left(r_{ij} \right),\vec{s}{j}\mathrm{and}\vec{F}{i} = \sum_{j}^{Neighb} \frac{\partial {J} \left(r_{ij} \right)}{ \partial r_{ij}} \left( \vec{s}{i}\cdot \vec{s}{j} \right) \vec{e}_{ij} $$

其中 $\hbar$ 是普朗克常数（metal 单位制下取约化形式，源码中由hbar = force->hplanck/MY_2PI计算，见 src/SPIN/pair_spin.cpp），$\vec{e}_{ij} = (\vec{r}_i - \vec{r}_j)/|\vec{r}_i-\vec{r}_j|$ 是两格点间的单位向量。当定义spin/exchange/biquadratic风格时，双二次项也会产生等价的力与磁力矩。

4.1 源码中的力学实现

compute_exchange_mech()实现了机械力计算 源码见。力的大小由对 $J(r_{ij})$ 求导得到：

Jex_mech = 1.0-ra-J2[itype][jtype]*ra*(2.0-ra); Jex_mech *= 8.0*Jex*rr*exp(-ra);

其中ra = rsq·iJ3、rr = sqrt(rsq)·iJ3，对应 $J'(r)$ 的闭式表达式；随后乘以自旋点积sdots = spi·spj与单位向量eij，得到沿径向的机械力分量。在compute()主循环中 源码见，每个自旋的磁力矩累加到atom->fm（磁力数组），机械力累加到atom->f（力学数组），f[i]加力、f[j]反向减力（在newton_pair on或j < nlocal时），能量通过ev_tally_xyz进行维里统计。这印证了该对势同时驱动自旋进动（经 fix nve/spin 辛积分）与晶格运动——这正是 Howto_spins 中所描述的"磁交换对自旋施加磁力矩，同时按自旋取向与原子间距对原子施力"的自旋-晶格耦合模型。

4.2 磁能输出

两种风格都维护一个emag数组用于逐原子磁能累积（emag[i] += evdwl），可通过compute spin输出总磁化、自旋温度与磁能等平均量 文档见，或通过compute property/atom输出每个原子的自旋分量与磁力分量 文档见。

五、pair_coeff 参数顺序与单位

对于spin/exchange风格，必须为每对原子类型通过pair_coeff命令（或由read_data/read_restart读取的 data 文件、restart 文件）按以下顺序定义系数：

对于spin/exchange/biquadratic风格，需要按以下顺序定义七个系数：

注意：$R_c$ 是该交换相互作用的截断半径，$a$、$b$、$d$ 是上文 $J(r_{ij})$ 函数的三个参数化系数（biquadratic风格中下标j、k分别对应 $J(r_{ij})$ 与 $K(r_{ij})$）。

所有系数均不可省略：如果未指定，spin/exchange风格将无法使用。源码中init_one()会检查setflag[i][j]，未设置全部 pair coeff 时直接报错"All pair coeffs are not set"源码见。此外，coeff()中参数读取逻辑为：arg[3]为 $R_c$，arg[4..6]为 $a,b,d$（biquadratic风格arg[7..9]为 $a_k,b_k,d_k$），arg[7]（或arg[10]）起可跟offset yes/no关键字，这也是文档示例中narg允许为 9（exchange）与 12（biquadratic）的原因。

六、offset 选项：消除磁能的机械残余压力

6.1 offset yes 时的能量替换

对于自旋-晶格模拟，对交换相互作用产生的机械力与能量进行偏移常常很有用，offset关键字即用于施加该偏移。当offset设为yes时，上面的能量定义替换为：

$$ H_{ex} = -\sum_{i,j}^N J_{ij} (r_{ij}) ,[ \vec{s}_i \cdot \vec{s}_j-1 ] $$

$$ H_{bi} = -\sum_{i, j}^{N} {J}{ij} \left(r{ij} \right), [ \vec{s}{i}\cdot \vec{s}{j} -1 ] -\sum_{i, j}^{N} {K}{ij} \left(r{ij} \right), [ \left(\vec{s}{i}\cdot \vec{s}{j}\right)^2 -1] $$

即把每个自旋对的自旋点积减去 1（或其平方减去 1），使完全平行排列（铁磁基态）的自旋对能量贡献为零。

6.2 关键性质与使用场景

该偏移只影响能量与机械力的计算，不修改进动向量的计算，因此不影响纯磁学性质。它保证了当所有自旋对齐时，磁能与相关的机械力（以及磁势产生的压力）为零。

官方文档特别提醒：在状态方程类计算（能量-体积、能量-压力曲线）中，该偏移项可能非常重要。设置offset可确保在晶体基态与平衡磁构型（典型为铁磁态）下压力为零；否则磁力可能产生残余压力，污染状态方程结果。当offset设为no时不做偏移，这也是默认选项（源码构造函数中e_offset = 0，见 src/SPIN/pair_spin_exchange.cpp）。

源码实现中，compute_energy()与compute_exchange_mech()都通过e_offset标志分支：e_offset == 1时用(sdots-1.0)（及(sdots*sdots-1.0)），e_offset == 0时用sdots（及sdots*sdots）源码见；非法值直接报错"Illegal option in pair exchange/biquadratic command"。

七、完整实战：构建可运行的自旋-晶格算例

仓库examples/SPIN提供了多个可运行算例。以下以 bcc 铁为例分析完整输入文件结构 in.spin.iron：

units metal atom_style spin dimension 3 boundary p p p atom_modify map array # 串行算法（sametag）所需 lattice bcc 2.8665 region box block 0.0 5.0 0.0 5.0 0.0 5.0 create_box 1 box create_atoms 1 box mass 1 55.845 set group all spin/random 31 2.2 # 随机初始化自旋，磁矩 2.2 μB velocity all create 100 4928459 rot yes dist gaussian pair_style hybrid/overlay eam/alloy spin/exchange 3.5 pair_coeff * * eam/alloy Fe_Mishin2006.eam.alloy Fe pair_coeff * * spin/exchange exchange 3.4 0.02726 0.2171 1.841 offset yes neighbor 0.1 bin neigh_modify every 10 check yes delay 20 fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0 # 外磁场（Zeeman） fix_modify 1 energy yes fix 2 all langevin/spin 0.0 0.0 21 # 自旋 Langevin 热浴 fix 3 all nve/spin lattice moving # 自旋-晶格辛积分 timestep 0.0001 compute out_mag all spin compute out_pe all pe compute out_ke all ke compute out_temp all temp variable magz equal c_out_mag[3] variable magnorm equal c_out_mag[4] variable emag equal c_out_mag[5] variable tmag equal c_out_mag[6] thermo_style custom step time v_magnorm v_tmag temp v_emag ke pe press etotal thermo 50 compute outsp all property/atom spx spy spz sp fmx fmy fmz dump 1 all custom 100 dump_iron.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3] run 1000

要点解读：

• atom_style spin是前提：SPIN 包全部pair/spin风格只有在构建时启用了 SPIN 包（参见 Build_package）且声明了atom_style spin时才可用。源码 pair_spin.cpp 在init_style()中强制检查：if (!atom->sp_flag) error->all(FLERR,"Pair spin requires atom/spin style")，并且要求newton pair on、使用完整邻居列表REQ_FULL，还会在未检测到fix nve/spin或fix neb/spin时给出警告。
• units metal强制：PairSpin::settings()校验if (strcmp(update->unit_style,"metal") != 0) error->all(... "Spin pair styles require metal units")源码见，spin 类对势只能在 metal 单位制下使用。
• hybrid/overlay组合：自旋交换项通常与 EAM 等原子间势叠加（铁用Fe_Mishin2006.eam.alloy，镍用Ni99.eam.alloy，钴用Co_PurjaPun_2012.eam.alloy），实现"原子势 + 磁性交换"的双通道力计算。
• offset yes的典型用法：铁算例在状态方程类分析中开启偏移，避免磁力产生残余压力。
• fix nve/spin lattice moving：moving表示自旋与晶格耦合运动；而 in.spin.bfo 中使用lattice frozen固定晶格，仅做纯自旋动力学——BFO 算例还同时叠加spin/magelec与spin/dmi三种自旋相互作用。
• 自旋初始化：set group all spin/random seed magmom以随机取向初始化自旋，也可用set group all spin mx my mz magmom指定方向。
• 磁量输出：compute spin的第 3/4/5/6 列分别对应 $M_z$、总磁化强度模、磁能与自旋温度，配合thermo_style custom实时监测；compute property/atom spx spy spz sp fmx fmy fmz输出逐原子自旋分量与磁力。

7.1 单元测试中的参考实现

仓库 unittest/force-styles/tests 提供了两种风格的数值参考：如 spin-pair-exchange.yaml 使用pair_style spin/exchange 3.5与pair_coeff * * exchange 3.4 0.02726 0.2171 1.841 offset yes，并记录 54 个原子的初末态力、磁力、应力与能量；spin-pair-exchange_biquadratic.yaml 则使用pair_style spin/exchange/biquadratic 4.0与pair_coeff * * biquadratic 4.0 0.05 0.03 1.48 0.05 0.03 1.48 offset yes。这些 YAML 测试文件可用于验证自编译 LAMMPS 在相同输入下的数值一致性（epsilon 级别为 5e-13）。

八、限制条件、相关命令与默认值

限制：所有pair/spin风格都属于 SPIN 包，只有在构建 LAMMPS 时启用该包（详见 Build_package）且声明atom_style spin时才可用；此外必须使用 metal 单位制、开启newton pair，并要求搭配fix nve/spin或fix neb/spin使用（源码 pair_spin.cpp 会发出警告）。

相关命令：atom_style spin、pair_coeff、pair_eam、以及配套的 fix nve/spin、min_spin 与自旋磁矩/磁能输出命令 compute spin。

默认值：offset关键字默认值为no（不施加偏移）。

参考文献：该实现与参数化方法基于 Tranchida, Plimpton, Thibaudeau 与 Thompson 在 Journal of Computational Physics, 372, 406-425 (2018) 发表的关于耦合磁自旋动力学与分子动力学的大规模并行辛算法论文，源码头部注释（pair_spin_exchange.cpp）亦要求引用该文献。

• 科研
• 科学计算
• 高性能计算