LAMMPS fix polarize 系列命令详解：DIELECTRIC 包中界面感应电荷的三种求解方法

科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载在含电解质界面的介电系统中如疏水壁、胶体、脂质双分子层界面上会自发产生感应电荷直接影响库仑相互作用与离子力。LAMMPS 的fix polarize/bem/gmres、fix polarize/bem/icc、fix polarize/functional三个 fix 命令隶属 DIELECTRIC 包就是专门用于在两个介电常数不同的介质界面处计算并更新感应电荷的工具。读完本文你将掌握这三个 fix 的语法与参数、atom_style dielectric数据文件的 15 列格式、fix_modify全部可调项、q_scaled输出机制并能基于仓库自带的 confined / nopbc 示例直接搭建可运行的界面极化模拟。一、命令语法与三种求解器三个 fix 共用同一套语法ID与group-ID的含义与fix命令一致fix ID group-ID style nevery tolerance参数说明ID、group-IDfix 的标识符与界面粒子所在原子组详见 fix 命令文档style取值为polarize/bem/gmres、polarize/bem/icc或polarize/functionalnevery该 fix 每隔多少个时间步被调用一次即ntimestep % nevery 0时才触发求解源码中pre_force有相应判断见 fix_polarize_bem_gmres.cpptolerance迭代求解器的相对收敛容差。文档给出的三个最简用法示例fix 2 interface polarize/bem/gmres 5 0.0001 fix 1 interface polarize/bem/icc 1 0.0001 fix 3 interface polarize/functional 1 0.0001三个求解器对应三种独立实现源码分别位于 fix_polarize_bem_gmres.cpp、fix_polarize_bem_icc.cpp 和 fix_polarize_functional.cppDIELECTRIC 包 README 对其定位如下polarize/bem/gmres基于边界元BEM方程、用 GMRES广义最小残差算法求解感应电荷方法出自 Barros 等人的论文polarize/bem/icc同样是边界元方程但用逐次超松弛successive over-relaxation / SOR迭代求解方法出自 Tyagi 等人的论文polarize/functional直接对能量泛函做变分优化方法出自 Jadhao 等人的论文。界面既可以是静止的如接触水的石墨烯层也可以是运动的如溶剂中的介电胶体——只要同时设置了合适的时间积分器例如fix rigid接口粒子会随运动更新fix 每个nevery步重新求解感应电荷。二、物理模型界面离散化与 dielectric 原子属性这三个 fix 计算的是两不可穿透介质界面处的感应电荷。界面必须离散化为一系列顶点vertices每个顶点对应一个边界元且被当作普通原子来存储和通信。因此必须使用atom_style dielectric见 atom_style 文档它额外定义了每个界面粒子的法向量、元素面积、局部介电失配等属性——从 atom_vec_dielectric.h 可以看到该原子风格新增的字段为mu法向量可通过mux/muy/muz像偶极矩分量一样访问、area、ed介电差、em介电均值、epsilon局部介电常数、curvature平均曲率和q_scaled。文档用一个水-二氧化硅体系说明了各属性的物理含义一侧介质介电常数为 78水另一侧为 4二氧化硅界面离散为 2000 个边界元每个界面试剂interface particle的法向量由二氧化硅指向水则法向介电差ed78 − 4 74平均介电值em(78 4) / 2 41每个边界元还带有自己的面积和局部平均曲率fix 用它计算局部电场中的修正项若要模拟带电界面界面试剂应带有非零电荷由其面积 × 表面电荷密度决定且其局部介电常数值取平均值em。对非界面粒子普通原子、带电粒子而言法向量、面积、介电差字段无意义且不被使用它们的epsilon字段会被内部使用计算库仑相互作用时把真实电荷除以该值进行缩放。例如模拟介电常数为 40 的隐式溶剂中携带 2 电荷的阳离子只需把q设为 2、epsilon设为 40无需手动缩放电荷任何带 dielectric 后缀的pair_style都会给出正确受力。文档同时指出一个隐含假设由于epsilon不随时间改变粒子在模拟中不能穿越界面。atom_style dielectric的 data 文件 Atoms 段包含15 列id mol type q x y z normx normy normz area/patch ed em epsilon curvature仓库示例 data.confined 中一个界面行如下1 0 1 0 0 0 9.99798 0 0 1 0.866 8 6 6 0而离子行data.confined 末尾为4001 0 2 1 15 20 15 0 0 1 1.0 8 6 10 0 4002 0 3 -1 25 20 25 0 0 1 1.0 8 6 10 0按 示例 README 的说明normx/normy/normz是该边界元界面处法单位向量的三个分量对真实离子无意义可填 0,0,1ed epsilon_out − epsilon_in法向从epsilon_in指向epsilon_outem (epsilon_out epsilon_in)/2界面粒子的epsilon取em离子的epsilon取所在介质的介电常数area/patch是元素面积对离子可填 1.0curvature为平均曲率球形界面取 1/半径平面取 0。三、感应电荷方程与三种算法3.1 边界元方程对polarize/bem/gmres与polarize/bem/icc界面组内顶点上的感应电荷密度的表达式为$$\sigma_b(\mathbf{s}) \frac{1 - \bar{\epsilon}}{\bar{\epsilon}}, \sigma_f(\mathbf{s}) - \epsilon_0 \frac{\Delta \epsilon}{\bar{\epsilon}}, \mathbf{E}(\mathbf{s})\cdot\mathbf{n}(\mathbf{s})$$其中$\sigma_b$界面试剂 $\mathbf{s}$ 处的感应电荷密度$\bar{\epsilon} (\epsilon_1 \epsilon_2)/2$界面试剂处平均介电常数即em$\Delta\epsilon \epsilon_1 - \epsilon_2$介电常数差即ed$\sigma_f$界面试剂处的自由电荷密度q/area$\mathbf{E}(\mathbf{s})$顶点处的电场$\mathbf{n}(\mathbf{s})$从 $\epsilon_2$ 介质指向 $\epsilon_1$ 介质的单位法向量。由于右端电场本身又依赖于感应电荷这是一个隐式方程需要迭代求解。在源码中右端向量 b 的构造见 compute_induced_charges()其中每个顶点计算double ndotE epsilon0e2q * (Ex * norm[i][0] Ey * norm[i][1] Ez * norm[i][2]) / epsilon[i]; double sigma_f q[i] / area[i]; buffer[idx] (1 - em[i]) * sigma_f - ed[i] * ndotE / MY_4PI;epsilon0e2q是单位换算因子把 $\epsilon_0 \cdot E$ 换算为电荷/面积²单位源码 setup() 中详细注释了real、metal、si、nano各单位制下该系数的推导并对不支持的单位制直接报错else if (strcmp(update-unit_style, lj) ! 0) error-all(FLERR, Only unit styles lj, real, metal, si and nano are supported);两个 BEM 求解器的迭代都在这两类终止条件之一满足时结束连续两次迭代间感应电荷的最大相对变化低于设定容差或迭代次数达到itr_max见下文fix_modify。polarize/bem/gmres采用 GMRES 算法Barros et al.在 gmres_solve() 中实现了 Arnoldi 正交化 Givens 旋转的标准流程。值得注意的是实现上并不显式存储矩阵 A而是把作用于向量封装为 apply_operator()把试探电荷写入q_scaled调用 pair/kspace 风格计算电场后读回残差则由 update_residual() 按 Barros 文献式 (60) 直接更新。收敛判据在 gmres_solve() 尾部rho tol_rel * normbBarros 文献建议的条件或rho rho_tol rho tol_abspolarize/bem/icc采用逐次超松弛算法Tyagi et al.。在 compute_induced_charges() 中每次迭代先用 pair/kspace 风格重新计算电场然后做 SOR 更新q_bound (1 - omega) * q_bound omega * ((1.0 / em[i] - 1) * q_free - (ed[i] / em[i]) * ndotE * area[i]);收敛判据为各顶点感应电荷相对变化的全局最大值MPI_Allreduce(..., MPI_MAX, ...)小于容差polarize/functional采用能量泛函变分法Jadhao et al.。从 源码头部注释 看它求解的是线性系统(Rww Rww^T) w q Rwq并自 2019 年 10 月起切换为共轭梯度conjugate gradient求解器。该 fix 会在init()中向 neighbor 请求完整邻居表NeighConst::REQ_FULLinit()其构造函数的容差是可选参数缺省为 1.0e-6构造函数。更完整的实现细节与推荐用法见 Nguyen 等人的论文(NguyenTD)Nguyen, Li, Bagchi, Solis, Olvera de la Cruz, Comput Phys Commun 241, 80-19 (2019)DIELECTRIC 包 README 也要求使用该包时引用该文献。3.2 依赖的 pair / kspace 风格文档要求同时使用带dielectric后缀的pair_style与kspace_style。源码中 fix 通过dynamic_cast逐一识别兼容风格update_efield() 列出的 pair 风格有lj/cut/coul/long/dielectric、lj/cut/coul/msm/dielectric、lj/cut/coul/cut/dielectric、lj/cut/coul/debye/dielectric、coul/long/dielectric、coul/cut/dielectric不匹配则报 Pair style not compatiblekspace 风格支持pppm/dielectric与msm/dielectric。电场计算按 pair 与 kspace 两部分分别累积到efield_pair、efield_kspace再求和。主要库仑风格文档见 pair lj/cut/coul/long/dielectric 与 kspace_style pppm/dielectric。四、q_scaled 输出与 restart 行为三个 fix 求出的感应电荷存储在q_scaled字段中。q字段始终是真实电荷、模拟中不改变对界面粒子q_scaled 真实电荷/局部介电常数epsilon 感应电荷对非界面粒子q_scaled 真实电荷/epsilon。GMRES 实现中这一点体现得很直接compute_induced_charges() 尾部if (induced_charge_idx[i] 0) { int idx induced_charge_idx[i]; q_scaled[i] induced_charges[idx] * area[i] q[i]; } else { q_scaled[i] q_backup[i]; }注意感应电荷induced_charges是电荷密度乘上area[i]才是该边界元的电荷量。文档给出的读取示例compute qs all property/atom q_scaled dump 1 all custom 1000 all.txt id type q x y z c_qs关于 restart该 fix 的任何信息都不写入二进制 restart 文件感应电荷会在setup()时重新求解初始化因此从 restart 恢复不影响正确性。polarize/bem/gmres与polarize/bem/icc还会计算一个 2 元素的全局向量global_freq 1size_vector 2构造函数可通过各类输出命令访问第一个元素是求解器终止时的迭代次数上限为itr_max第二个元素是 RMS 误差。compute_vector()的实现见 fix_polarize_bem_gmres.cpp#L969-L977可用variable x1 c_3[1]/c_3[2]3为 fix ID配合thermo_style custom或fix ave/time监控收敛质量。polarize/functional则不提供向量输出。五、fix_modify 参数详解通过 fix_modify 命令 可修改以下设置三个 fix 的modify_param实现分别见 GMRES 与 ICCfix_modify ID itr_max arg fix_modify ID dielectrics ediff emean epsilon area charge fix_modify ID kspace arg fix_modify ID rand max seed fix_modify ID mr arg fix_modify ID omega argitr_max argarg 每步求解的最大迭代次数控制收敛迭代的上界dielectrics ediff emean epsilon area charge覆盖group-ID中所有原子的介电属性。ediff 介电差emean 介电均值epsilon 局部介电值area 元素面积charge 界面真实电荷任一参数传NULL则保留 data 文件中的现有值。注意该设置作用于整个组——若两个不同界面需要不同属性需要为每个界面组分别建立一对fixfix_modify。源码中NULL判断与正数校验逻辑见 set_dielectric_params()kspace argyes或no开关长程相互作用参与感应电荷计算默认 yes即若存在 kspace 风格则自动启用若指定 yes 但系统没有 kspace 风格会打印警告 No Kspace pair style available 并禁用见 setup()rand max seed以均匀分布生成范围 −max/2 到max/2 的随机初始感应电荷随机数种子为seed用于测试收敛性。源码中实际会进一步把随机电荷全局减均值使总电荷为零init()mr arg仅polarize/bem/gmres求解表面电荷时使用的 q-向量Krylov 空间最大数量omega arg仅polarize/bem/iccTyagi 文献定义的松弛参数一般应设在 0 到 2 之间。默认值汇总参数默认值适用itr_max50全部kspaceyes全部omega0.7仅 ICCmr组内原子数 − 1仅 GMRES随机电荷初始化不启用全部源码印证itr_max 50在 GMRES 构造函数、ICC 构造函数 中初始化mr默认在init()中被钳制为mat_dim - 1fix_polarize_bem_gmres.cpp#L202-L203omega 0.7见 fix_polarize_bem_icc.cpp#L73。另外文档提醒局部介电常数epsilon也可以独立用 set 命令 设置示例脚本中正是这么做的。六、完整示例examples/PACKAGES/dielectric文档指引的示例位于examples/PACKAGES/dielectric目录包含两套场景见 示例 READMEdata.confined两个等量异号点电荷被限制在两个平行界面之间epsilon12 | epsilon210 | epsilon12data.sphere两个点电荷位于球形界面试剂之外epsilon_in1epsilon_out10。6.1 confined 示例周期性边界 PPPM 长程入口脚本 in.confined.gmres 只有两行通过变量选择求解器后 include 公共脚本variable method index gmres # gmres BEM/GMRES # icc BEM/ICC* # dof Direct optimization of the functional include confined.incin.confined.icc 与 in.confined.dof 只需把method改为icc或dof。核心逻辑在 confined.inc完整流程如下# 界面法向量方向须与 ed 一致从 epsilon1 区指向 epsilon2 区 variable epsilon1 index 20 variable epsilon2 index 10 variable data index data.confined newton off units lj atom_style dielectric atom_modify map array dimension 3 boundary p p f # 根据介电常数在脚本内计算界面试剂属性 variable ed equal v_epsilon2 - v_epsilon1 variable em equal (v_epsilon2 v_epsilon1)/2 variable epsilon equal 1.0 # 当前未使用的 patch epsilon variable area equal 0.866 # patch 面积与 data 文件一致 read_data ${data} group interface type 1 group ions type 2 3 group cations type 2 group anions type 3 # 设置离子所在介质的介电常数epsilon 用 set 命令独立设置 set group cations epsilon ${epsilon2} set group anions epsilon ${epsilon2} pair_style lj/cut/coul/long/dielectric 1.122 10.0 pair_coeff * * 1.0 1.0 pair_coeff 1 1 0.0 1.0 kspace_style pppm/dielectric 0.0001 kspace_modify slab 3.0 neigh_modify every 1 delay 0 check yes one 5000 fix 1 ions nve # 用 fix / fix_modify 设置界面组属性并按 method 选择求解器 if ${method} gmres then fix 3 interface polarize/bem/gmres 1 1.0e-4 fix_modify 3 itr_max 50 dielectrics ${ed} ${em} ${epsilon} ${area} NULL elif ${method} icc fix 3 interface polarize/bem/icc 1 1.0e-4 fix_modify 3 itr_max 50 dielectrics ${ed} ${em} ${epsilon} ${area} NULL elif ${method} dof fix 3 interface polarize/functional 1 0.0001 fix_modify 3 dielectrics ${ed} ${em} ${epsilon} ${area} NULL else print Unsupported polarization solver thermo 1000 thermo_style custom step evdwl ecoul elong epair #f_3 thermo_modify flush yes run 0该脚本展示了文档中提到的全部关键要素newton offfix 内部假设 pair 使用完整邻居表并自行做通信、atom_modify map array配合 fix 内按 atom tag 建立tag2mat映射、set命令设置离子epsilon、fix_modify dielectrics用变量覆盖界面属性。confined.inc 中还保留了被注释掉的调试输出可用compute efield/atom文档导出离子与界面粒子上的电场分量来验证极化分布#compute ef all efield/atom #dump 1 all custom 100 all.dump id mol type q x y z #fx fy fz c_ef[1] c_ef[2] c_ef[3]6.2 nopbc 示例非周期边界 大截断nopbc.inc 演示无长程求解器的路线boundary f f fpair_style lj/cut/coul/cut/dielectric 1.122 20.0使用足够大的截断来模拟长程neigh_modify one 5000保证大截断下邻居表完整。三种求解器的选择方式与 confined 相同。README 给出的运行方式为需已编译 DIELECTRIC 包mpirun -np 8 ../../../src/lmp_mpi -in in.confined仓库中同时保留着三种求解器分别的运行日志如log.*.confined.gmres.g.4可直接对比同一场景下不同算法的收敛迭代数与能量。6.3 单元测试除示例外unittest/force-styles/tests/in.dielectric 配合 data.dielectric 被各 dielectric pair/kspace 风格的 YAML 测试如 kspace-pppm_dielectric.yaml、mol-pair-lj_cut_coul_long_dielectric.yaml复用是验证atom_style dielectric与 dielectric 风格受力一致性的回归基准修改介电相关代码时可参考。七、构建前提与使用限制7.1 包依赖这些 fix 属于DIELECTRIC 包只有 LAMMPS 用该包编译后才可用且必须同时安装 KSPACE 包构建细节见 Build package 文档。依赖关系在 src/DIELECTRIC/Install.sh 中被硬校验# compute efield/atom is only usable when all styles # are installed, which in turn requires KSPACE if (test $1 1) then if (test ! -e ../pppm.cpp) then echo Must install KSPACE package with DIELECTRIC package exit 1 fi if (test ! -e ../pair_lj_cut_coul_debye.cpp) then echo Must install EXTRA-PAIR package with DIELECTRIC package exit 1 fi fi即安装 DIELECTRIC 前系统里必须已有 KSPACEpppm.cpp与 EXTRA-PAIRpair_lj_cut_coul_debye.cpp文件。7.2 单位制与功能限制polarize/bem/gmres与polarize/bem/icc目前只支持lj、real、metal、si、nano五种 units源码对其它单位制直接报错见上文 setup() 代码polarize/functional尚不支持带电界面界面粒子q非零的场景模型假设粒子不穿越界面因为epsilon属性在模拟中不改变两个 BEM fix 的fix_modify属性覆盖是组级别的同一 fix 管理的界面共享同一组ediff/emean/epsilon/area/charge。7.3 相关命令与参考资料相关命令pair_coeff、fix polarize、read_data、pair_style lj/cut/coul/long/dielectric、kspace_style pppm/dielectric、compute efield/atom输出机制见 Howto output 与 restart。方法学参考文献与文档对应(Barros)Barros, Sinkovits, Luijten, J. Chem. Phys, 140, 064903 (2014)——BEM 方程与 GMRES 求解(Tyagi)Tyagi, Suzen, Sega, Barbosa, Kantorovich, Holm, J Chem Phys, 132, 154112 (2010)——逐次超松弛ICC算法及omega参数(Jadhao)Jadhao, Solis, Olvera de la Cruz, J Chem Phys, 138, 054119 (2013)——能量泛函变分法(NguyenTD)Nguyen, Li, Bagchi, Solis, Olvera de la Cruz, Comput Phys Commun 241, 80-91 (2019)——本实现的完整描述与推荐用法使用 DIELECTRIC 包发表论文时应引用。赞分享科研科学计算高性能计算【免费下载链接】lammpsPublic development project of the LAMMPS MD software package项目地址https://gitcode.com/gh_mirrors/la/lammps点击查看免费下载相关推荐LAMMPS 中 compute efield/atom 命令详解使用 DIELECTRIC 包计算原子级电场LAMMPS 中 compute efield/atom 命令详解使用 DIELECTRIC 包计算原子级电场 本文以 LAMMPS 官方文档 doc/src科研科学计算高性能计算LAMMPS compute dipole 命令详解:体系偶极矩计算与 TIP4P 电荷位点修正LAMMPS compute dipole 命令详解:体系偶极矩计算与 TIP4P 电荷位点修正 本文围绕 LAMMPS 的 compute dipole 与科研科学计算高性能计算LAMMPS 介电常数设置指南dielectric 命令的语法、物理含义与源码实现LAMMPS 介电常数设置指南dielectric 命令的语法、物理含义与源码实现 dielectric 是 LAMMPS 中用于全局调整库仑Coulomb科研科学计算高性能计算上一篇GolangBypassAV 项目使用教程下一篇SignTools 开源项目教程创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考