【LAMMPS学习】八、基础知识(5.8)LAMMPS 中热化 Drude 振荡器教程

8. 基础知识

此部分描述了如何使用 LAMMPS 为用户和开发人员执行各种任务。术语表页面还列出了 MD 术语，以及相应 LAMMPS 手册页的链接。 LAMMPS 源代码分发的 examples 目录中包含的示例输入脚本以及示例脚本页面上突出显示的示例输入脚本还展示了如何设置和运行各种模拟。

8.1.通用基础知识

8.2. 设置入门

8.3. 分析入门

8.4. 力场入门

8.5. 软件包入门

8.5.1.有限尺寸球形和非球形粒子

8.5.2. 粒度模型

8.5.3.体粒子

8.5.4.黏合粒子(BMP)模型

8.5.5.极化模型

8.5.6. 绝热核/壳模型

8.5.7.Drude感应偶极子

8.5.8. LAMMPS 中热化 Drude 振荡器教程

本教程介绍如何在 LAMMPS 中使用 Drude 振荡器，通过 DRUDE 包模拟可极化系统。作为说明，记录了 250 个苯酚分子模拟的输入文件。首先，必须在激活 DRUDE 包的情况下编译 LAMMPS。然后，必须修改数据文件和输入脚本以包含德鲁德偶极子以及如何处理它们。

可用的示例输入脚本：examples/PACKAGES/drude

Drude感应偶极子概述

可极化原子在外部电场（例如由周围粒子产生的电场）的作用下获得感应电偶极矩。德鲁德振荡器通过两个固定电荷表示这些偶极子：核心 (DC) 和由谐波势束缚的德鲁德粒子 (DP)。德鲁德粒子可以被认为是电子云，其中心可以偏离相应原子核的位置。

核心-德鲁德对的质量总和应该是初始（未分裂）原子 𝑚𝐶+𝑚𝐷=𝑚 的质量。它们的电荷总和应该是初始（未分裂）原子的电荷 𝑞𝐶+𝑞𝐷=𝑞 。核心和 Drude 伙伴之间应存在谐波势，力常数 𝑘𝐷 且平衡距离为零。 harmonic bond 𝐾𝐷=𝑘𝐷/2 的（半）刚度和德鲁德电荷 𝑞𝐷 与原子极化率 𝛼 的关系为：

$K_D=\frac{1}{2} \frac{q_D^2}{\alpha}$

理想情况下，德鲁德粒子的质量应该很小，调和键的刚度应该很大，以便德鲁德粒子保持靠近核心。可以按照不同的策略选择 Drude 质量、Drude 电荷和力常数的值，如以下极化力场示例所示：

Lamoureux 和 Roux 建议采用全局半刚度， 𝐾𝐷 = 500 kcal/(mol Ang 2 ) - 对应于力常数 𝑘𝐷 = 4184 kJ /(mol Ang 2 ) - 对于所有类型的核心-Drude 键，对于所有类型的 Drude 粒子，全局质量 𝑚𝐷 = 0.4 g/mol（或 u），并且使用方程 (1) 根据原子极化率计算各个原子类型的德鲁德电荷。极化 CHARMM 力场遵循这一选择。
另外， Schroeder 和 Steinhauser 建议对所有 Drude 粒子采用全局电荷 𝑞𝐷 = -1.0e 和全局质量 𝑚𝐷 = 0.1 g/mol（或 u），并计算力方程（1）中每种类型的核心-德鲁德键的常数。这些作者使用的时间步长在 0.5 到 2 fs 之间，Drude 振荡器的自由度保持在 1 K。
在这两个力场中，氢原子都被视为不可极化的。

德鲁德粒子的运动可以通过迭代、自洽的过程最小化每个时间步长的感应偶极子的能量来计算。德鲁德粒子可以是无质量的，因此不会贡献动能。然而，宽松的方法计算速度慢。扩展拉格朗日方法可用于计算德鲁德粒子的位置，但这要求它们具有质量。在这种情况下，将与德鲁德振子相关的自由度与普通原子的自由度解耦非常重要。在与模拟的其余部分相当的温度下热化 Drude 偶极子会导致几个问题（动能转移、非常短的时间步长等），这些问题可以通过“冷 Drude”技术(Lamoureux 和 Roux)来解决。

两个密切相关的模型用于通过“弹簧上的电荷”来表示极化：核壳模型和德鲁德模型。尽管基本思想相同，但核壳模型通常用于离子/晶体材料，而德鲁德模型通常用于分子系统和流体状态。在离子晶体中，每个离子周围的对称性以及它们之间的距离使得核-壳模型足够稳定。但为了适用于分子/共价系统，Drude 模型包括两个重要特征：

根据德鲁德粒子相对于其核心的简化坐标，可以在非常低的温度下恒温与感应偶极子相关的附加自由度。这使得轨迹接近松弛感应偶极子的轨迹。
由于距离短，共价键原子上的德鲁德偶极子相互作用太强，因此原子可能会捕获邻居的德鲁德粒子（壳），或者同一分子内的诱导偶极子可能会对齐太多。为了避免这种情况，可以通过 Thole函数来对构成感应偶极子的点电荷之间的相互作用进行阻尼。

数据文件的准备

该数据文件类似于atom_style full 的标准LAMMPS 数据文件。 DP 以及将它们连接到 DC 的谐波键应作为正常原子和键出现在数据文件中。

您可以使用极化器工具（随 DRUDE 包分发的 Python 脚本）将非极化数据文件（此处为 data.102494.lmp）转换为极化数据文件（data-p.lmp）

polarizer -q -f phenol.dff data.102494.lmp data-p.lmp

这将自动插入新的原子和键。 DC 和 DP 的质量和电荷是根据 pheno.dff 以及 DC-DP 键常数计算的。文件 pheno.dff 包含原子类型的极化率和 Drude 粒子的质量，例如：

# units: kJ/mol, A, deg
# kforce is in the form k/2 r_D^2
# type  m_D/u   q_D/e    k_D   alpha/A3  thole
OH      0.4    -1.0    4184.0   0.63     0.67
CA      0.4    -1.0    4184.0   1.36     2.51
CAI     0.4    -1.0    4184.0   1.09     2.51

该文件中不存在氢原子，因此它们将被视为不可极化原子。在非极化数据文件 data.102494.lmp 中，必须将与原子类型编号相对应的原子名称指定为 Masses 部分行末尾的注释。您可能需要编辑它才能添加这些名称。它应该看起来像

Masses1 12.011 # CAI
2 12.011 # CA
3 15.999 # OH
4 1.008  # HA
5 1.008  # HO

基本输入文件

原子样式应设置为（或派生自）完整，以便您可以定义原子电荷和分子键、角度、二面角......

偏振器工具还输出与输入脚本相关的某些行（这些行的使用将在下面解释）。为了让 LAMMPS 识别出您正在使用 Drude 振荡器，您应该使用修复 drude。命令是

fix DRUDE all drude C C C N N D D D

drude 关键字后面的 N、C、D 具有以下含义：每种原子类型都有一个标记。该标签是 C（代表 DC）、D（代表 DP）和 N（代表不可极化原子）。这里，原子类型 1 至 3（C 和 O 原子）是 DC，原子类型 4 和 5（H 原子）是不可极化的，原子类型 6 至 8 是新创建的 DP。

通过识别修复德鲁德，LAMMPS 将找到并存储匹配的 DC-DP 对，并将 DP 视为与特殊键关系中的 DC 等效。可能需要扩展存储这种特殊关系的空间。在这种情况下，应使用 read_data 和create_box 命令的 extra/special/per/atom 关键字来保留额外的空间。对于我们的苯酚，还有 1 个特殊的邻居需要空间。否则 LAMMPS 会崩溃并给出所需的值。

read_data data-p.lmp extra/special/per/atom 1

假设我们要在 300 K 下运行简单的 NVT 仿真。请注意，Drude 振荡器需要在低温下热化才能近似自洽场 (SCF)，因此不可能仿真 NVE 系综有了这个包。由于偶极子近似为带电 DC-DP 对，因此pair_style 必须包括库仑相互作用，例如 lj/cut/coul/long 与 kspace_style pppm。例如，截断值为 10.，精度为 1.e-4：

pair_style lj/cut/coul/long 10.0
kspace_style pppm 1.0e-4

现在对于热化，最简单的选择是使用fix langevin/drude。

这将温度为 300 的 Langevin 恒温器应用于 DC-DP 对的质心，弛豫时间为 100，随机种子为 12345。此修复也将温度为 1 的 Langevin 恒温器应用于 DP 周围的相对运动它们的 DC，弛豫时间为 20，随机种子为 13977。只有 DC 和不可极化原子需要位于此修复组中。 LAMMPS 将对 DP 及其 DC 进行恒温。为此，幽灵原子需要知道它们的速度。因此您需要添加以下命令：

comm_modify vel yes

为了避免由于 Langevin 恒温器对 DC 的随机力而导致整个系统的质心发生漂移，您可以在修复线末尾添加零是选项。

如果使用修复抖动来约束 C-H 键，则应在修复 langevin/drude 之后调用它，以获得更高的准确性。

fix SHAKE ATOMS shake 0.0001 20 0 t 4 5

注： fix shake的组不得包括 DP。如果组 ATOMS 是由非 DPs 原子类型定义的，您可以使用

由于fix langevin/drude不执行时间积分（仅修改力，但不更新位置/速度），因此应结合使用修复 nve。

fix NVE all nve

为了避免原子逐渐冻结并且整个系统的质心漂移得越来越快的飞行冰块伪影，可以使用固定动量。例如：

fix MOMENTUM all momentum 100 linear 1 1 1

最后，如果您在 dump_modify … element … 命令中使用原子类型元素，请不要忘记通过添加 DP 的元素类型来更新它们。这里例如

dump DUMP all custom 10 dump.lammpstrj id mol type element x y z ix iy iz
dump_modify DUMP element C C O H H D D D

输入文件现在应该可以使用了！
您会注意到 LAMMPS 计算的全局温度 Thermo_temp 不是想要的 300.K。这是因为 LAMMPS 在默认计算中将 DP 视为标准原子。如果要输出 DC-DP 对质心和 DP 相对于其 DC 的温度，则应使用compute temp_drude

compute TDRUDE all temp/drude

然后使用 Thermo_style 自定义分别使用 c_TDRUDE[1] 和 c_TDRUDE[2] 输出 Drude 振荡器的正确温度。这些值平均应接近 300.0 和 1.0。

thermo_style custom step temp c_TDRUDE[1] c_TDRUDE[2]

Thole 筛选

弹簧上的点电荷所代表的偶极相互作用可能不稳定，例如，如果原子极化率太高，则DP可能会逃离其DC并被另一个DC捕获，这使得力和能量发散并导致模拟崩溃。即使没有达到这种极端情况，同一分子上附近偶极子之间的相关性也可能被夸大。通常，特殊的键关系会阻止键合的相邻原子看到彼此 DP 的电荷，因此问题并不总是出现。通过使用 *pair_style thole*可以使用屏蔽偶极子-偶极子相互作用。这是作为对库仑对样式的校正来实现的，它在短距离内抑制了代表感应偶极子的电荷之间的相互作用。它可与任何标准涂层对样式混合/覆盖。在我们的示例中，我们将使用

pair_style hybrid/overlay lj/cut/coul/long 10.0 thole 2.6 10.0

这告诉 LAMMPS 我们正在使用两个pair_style。第一个如上（lj/cut/coul/long 10.0）。第二个是一个整体的pair_style，默认筛选因子为2.6（Noskov），截止值为10.0。

由于混合/叠加不支持混合规则，因此应显式定义 i < j 的所有原子类型对的相互作用系数。偏振器脚本的输出可用于完成输入文件的pair_coeff 部分。在我们的示例中，这将如下所示：

pair_coeff    1    1 lj/cut/coul/long    0.0700   3.550
pair_coeff    1    2 lj/cut/coul/long    0.0700   3.550
pair_coeff    1    3 lj/cut/coul/long    0.1091   3.310
pair_coeff    1    4 lj/cut/coul/long    0.0458   2.985
pair_coeff    2    2 lj/cut/coul/long    0.0700   3.550
pair_coeff    2    3 lj/cut/coul/long    0.1091   3.310
pair_coeff    2    4 lj/cut/coul/long    0.0458   2.985
pair_coeff    3    3 lj/cut/coul/long    0.1700   3.070
pair_coeff    3    4 lj/cut/coul/long    0.0714   2.745
pair_coeff    4    4 lj/cut/coul/long    0.0300   2.420
pair_coeff    *    5 lj/cut/coul/long    0.0000   0.000
pair_coeff    *   6* lj/cut/coul/long    0.0000   0.000
pair_coeff    1    1 thole   1.090   2.510
pair_coeff    1    2 thole   1.218   2.510
pair_coeff    1    3 thole   0.829   1.590
pair_coeff    1    6 thole   1.090   2.510
pair_coeff    1    7 thole   1.218   2.510
pair_coeff    1    8 thole   0.829   1.590
pair_coeff    2    2 thole   1.360   2.510
pair_coeff    2    3 thole   0.926   1.590
pair_coeff    2    6 thole   1.218   2.510
pair_coeff    2    7 thole   1.360   2.510
pair_coeff    2    8 thole   0.926   1.590
pair_coeff    3    3 thole   0.630   0.670
pair_coeff    3    6 thole   0.829   1.590
pair_coeff    3    7 thole   0.926   1.590
pair_coeff    3    8 thole   0.630   0.670
pair_coeff    6    6 thole   1.090   2.510
pair_coeff    6    7 thole   1.218   2.510
pair_coeff    6    8 thole   0.829   1.590
pair_coeff    7    7 thole   1.360   2.510
pair_coeff    7    8 thole   0.926   1.590
pair_coeff    8    8 thole   0.630   0.670

对于thole对类型，系数为

以立方长度为单位的原子极化率
Thole函数的筛选因子（可选，默认值由pair_style命令指定）
截止值（可选，默认值由pair_style命令定义）

特殊邻居具有由special_bonds 命令的库因子（上例中的0.0、0.0 和0.5）筛选的电荷-电荷和电荷-偶极子相互作用。在不使用pair_style thole的情况下，偶极子-偶极子相互作用通过相同的因子进行筛选。通过使用pair_style thole，偶极子-偶极子相互作用可以通过Thole函数进行筛选，无论它们的特殊关系如何（当然每个DC-DP对内除外）。例如，考虑 1-2 个邻居：使用pair_style thole，它们的偶极子将看到彼此（尽管 coul 因子为 0），并且这些偶极子之间的相互作用将被 Thole 函数抑制。

恒温器和恒压器

将 Nose-Hoover 恒压器与 langevin/drude 恒温器结合使用非常简单，使用 fix nph 而不是 nve。例如：

fix NPH all nph iso 1. 1. 500

也可以使用 Nose-Hoover 代替 Langevin 恒温器。这需要在时间积分fixes之前和之后使用 *fix drude/transform* 。修复 drude/transform/direct 将原子质量、位置、速度和力转换为简化表示，其中 DC 转换为 DC-DP 对的质心，DP 转换为其相对于 DC 的相对位置。fix drude/transform/inverse执行逆向变换。对于 NVT 模拟，DC 和原子相对于 DC 的温度为 300 K，DP 的温度为 1 K，则可以使用

fix DIRECT all drude/transform/direct
fix NVT1 ATOMS nvt temp 300. 300. 100
fix NVT2 DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse

对于我们的苯酚示例，这些分组将定义为

group ATOMS  type 1 2 3 4 5 # DCs and non-polarizable atoms
group CORES  type 1 2 3     # DCs
group DRUDES type 6 7 8     # DPs

请注意，使用 fixes drude/transform 时，不需要指定 comm_modify vel yes，因为fixes无论如何都会这样做（多次并且对于力）。

使用 Nose-Hoover 恒压器和恒温器运行 NPT 模拟有点棘手。首先，体应该只积分一次。因此，DC 和原子的fixs应该是 npt，而 DP 的fixs应该是 nvt（反之亦然）。其次，fixs npt 计算全局压力，从而计算全局温度，无论fixs组如何。我们确实希望压力对应于整个系统，但我们希望温度仅对应于固定组。然后我们必须为此使用 fix_modify 命令。最后，恒温和恒压指令块将如下所示

compute TATOMS ATOMS temp
fix DIRECT all drude/transform/direct
fix NPT ATOMS npt temp 300. 300. 100 iso 1. 1. 500
fix_modify NPT temp TATOMS press thermo_press
fix NVT DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse

热化 Drude 模型的另一种选择是使用 (Son) 提出的温度分组 Nose-Hoover (TGNH) 恒温器。这是通过 fix tgnvt/drude 和 fix tgnpt/drude 来实现的。它将动能分为三个部分：分子质心 (COM) 运动、原子或原子-德鲁德对相对于分子 COM 的运动以及原子-德鲁德对的相对运动。独立的 Nose-Hoover 链适用于每种类型的运动。当使用TGNH时，分子、原子和Drude运动的温度可以用 thermo_style 命令打印出来。

使用 TGNH 恒温器进行 NVT 模拟

comm_modify vel yes
fix TGNVT all tgnvt/drude temp 300. 300. 100 1. 20
thermo_style custom f_TGNVT[1] f_TGNVT[2] f_TGNVT[3]

使用 TGNH 恒温器进行 NPT 模拟

comm_modify vel yes
fix TGNPT all tgnpt/drude temp 300. 300. 100 1. 20 iso 1. 1. 500
thermo_style custom f_TGNPT[1] f_TGNPT[2] f_TGNPT[3]

刚体

您可能想将分子模拟为刚体（但可极化）。常见的例子有SWM4-NDP等水模型，它是一种可极化的TIP4P水。刚体和 DP 应该单独集成，即使使用 Langevin 恒温器也是如此。让我们回顾一下不同的恒温器和整体组合。

使用 Langevin 恒温器的 NVT 系综：

comm_modify vel yes
fix LANG all langevin/drude 300. 100 12435 1. 20 13977
fix RIGID ATOMS rigid/nve/small molecule
fix NVE DRUDES nve

使用 Nose-Hoover 恒温器的 NVT 系综：

fix DIRECT all drude/transform/direct
fix RIGID ATOMS rigid/nvt/small molecule temp 300. 300. 100
fix NVT DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse

带 Langevin 恒温器的 NPT 系综：

comm_modify vel yes
fix LANG all langevin/drude 300. 100 12435 1. 20 13977
fix RIGID ATOMS rigid/nph/small molecule iso 1. 1. 500
fix NVE DRUDES nve

使用 Nose-Hoover 恒温器的 NPT 系综：

compute TATOM ATOMS temp
fix DIRECT all drude/transform/direct
fix RIGID ATOMS rigid/npt/small molecule temp 300. 300. 100 iso 1. 1. 500
fix_modify RIGID temp TATOM press thermo_press
fix NVT DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse

参考文献

(Lamoureux and Roux) Lamoureux and Roux, J Chem Phys, 119, 3025-3039 (2003)

(Schroeder) Schroeder and Steinhauser, J Chem Phys, 133, 154511 (2010).

(Thole) Chem Phys, 59, 341 (1981).

(Noskov) Noskov, Lamoureux and Roux, J Phys Chem B, 109, 6705 (2005).

(SWM4-NDP) Lamoureux, Harder, Vorobyov, Roux, MacKerell, Chem Phys Let, 418, 245-249 (2006)

(Son) Son, McDaniel, Cui and Yethiraj, J Phys Chem Lett, 10, 7523 (2019).