如何在VASP中约束磁矩¶
本文是 Tharindu 的一篇讲述如何在 VASP 中约束磁矩的文章。Tharindu 在测试的时候用到了 VASP 6.4.0 以上的版本，但是对于我们使用其他 VASP 版本进行约束磁矩计算也很有启发。
原文链接：How to constrain magnetic moments in VASP
一些从头计算需要约束离子磁矩。例如，在确定各向异性的方向（平面内、平面外等）或进行 4 态映射以计算磁各向异性（各向异性 J、DMI、单离子各向异性等）时。在本教程中，我将总结我从 VASP 文档中学到的内容以及我的个人经验，希望能让你在自己的计算中有一个良好的开端。
一个简单的起点是使用 MAGMOM 标签来定义每个离子的初始磁矩。然而，最终系统的离子可能会偏离这种初始配置以获得最低的能量。为了解决这个问题，可以在非共线计算中限制磁矩。回想一下，非共线计算使用 vasp_ncl 可执行文件而不是 vasp_std，并使用标签 LNONCOLLINEAR = .TRUE. 和 LSORBIT = .TRUE.（打开自旋轨道耦合），而不是 ISPIN=2（共线自旋计算）。据我所知，如果你想在共线自旋计算中约束力矩（你只关心向上自旋 +1 或向下自旋 -1，而不是一般的向量），你仍然需要通过非共线设置。
对于非共线计算，每个离子的磁矩应指定为与笛卡尔坐标系中的向量相对应的三元数组。例如，如果我有一个含有两个离子的系统，我的第一个离子需要（1,2,0），第二个需要（3,1,1），磁矩应该设置为：
这也是M_CONSTR 标签所需的相同格式（稍后讨论）。如果你正在模拟具有向上、向下、向下自旋的 3 原子系统的共线自旋计算，如果您在 2D 材料中具有平面外各向异性，则将需要使用：
如果你想快速确定各向异性的方向（平面内与平面外等），一个简单的测试是检查各种配置。
对于 x 面内：
对于 y 面内：
对于 z 面内：
使用矢量 \(\displaystyle\frac{1}{\sqrt{2}}(1,1,0)\) 在平面内 xy 平面内与 x 轴成 45 度角：
要实际约束磁矩，需要在设置非共线计算标记后使用 I_CONSTRAINED_M 标记。
I_CONSTRAINED_M=1 仅约束方向，I_CONSTRAINED_M=2 约束大小和方向，I_CONSTRAINED_M=4 约束方向和符号。
注：至少需要 VASP 6.4.0 版本才能使用 I_CONSTRAINED_M=4。
我可能错了，但根据我的经验，对于 I_CONSTRAINED_M=1，VASP 似乎无法区分沿同一平行线的方向。例如，无法区分 0 0 1 和 0 0 -1 。但是 I_CONSTRAINED_M=4 区分了这两者，并将沿着包括符号在内的给定方向进行精确约束。
例如，我使用 \(\rm Nb_3Cl_8Nb\)，其中 \(\rm Nb\) 离子被视为磁性，而 \(\rm Cl\) 离子则不是。
I_CONSTRAINED_M = 4
RWIGS = 1.270 1.111
LAMBDA = 9
MAGMOM = 1 0 0 1 0 0 1 0 0 24*0
M_CONSTR = 1 0 0 1 0 0 1 0 0 24*0
在 INCAR 中，对于 \(\rm Cl\) 离子，我写了 24*0 而不是写了 8 次 0 0 0。
还要注意，我将 MAGMOM 和 M_CONSTR 设置为相同。M_CONSTR 指定所需的时刻，而 MAGMOM 给出初始时刻。根据我的经验，我们需要将两个标签设置为相同，这对我们获得所需的结果非常重要。
此外，请注意参数 RWIGS 和 LAMBDA。这些标签的一个恼人的地方是，必须进行系统的收敛测试来确定它们的最佳值。
对于 RWIGS，Wigner-Seitz 半径（单位：\(\rm \AA\)），VASP 文档给出了一些确定最佳值的建议（对于非单原子系统，无法明确定义）。然而，这些建议可能很耗时，因此我有时会使用 POTCAR 文件中给出的 RWIGS 值。在这种情况，命令行运行 grep -e RWIGS POTCAR，得到：
RWIGS = 2.400; RWIGS = 1.270 wigner-seitz radius (au A)
RWIGS = 2.810; RWIGS = 1.487 wigner-seitz radius (au A)
根据 RWIGS 估计 \(\rm \AA\) 值，我分别为 \(\rm Nb\) 和 \(\rm Cl\) 选择了 1.270 和 1.487。另一种方法（在文档中列出）是使用像这样的数据集1中给出的原子共价半径。但是，请注意，不同的来源可能会对此给出不同的值。最终重要的是，计算出的结果是有意义的。
作为一种工作流程策略，我将首先使用 POTCAR 中的 RWIGS 值，并使用这些值进行 LAMBDA 参数的测试。然后，我将使用 LAMBDA 的最佳值来测试 RWIGS 的其他值，并从进行微调，同时确保积分磁矩与预期一致（见下一段）。
对于 LAMBDA，此参数对应于添加到系统中的能量损失。为了沿特定方向约束磁矩，必须添加能量以迫使磁矩约束。I_CONSTRAINED_M 手册规定，应通过逐渐增加 LAMBDA 来获得总 DFT 能量，直到惩戒函数对能量做的贡献变得非常小。
您可以在 OSZICAR 文件末尾的最后一个块中找到此贡献的值（我推测其单位为 eV），其中它们给出：
E_p = 0.36856E-07 lambda = 0.500E+02
<lVp>= 0.30680E-02
DBL = -0.30680E-02
ion MW_int M_int
1 -0.565 0.000 0.000 -0.770 0.000 0.000
2 0.565 0.000 0.000 0.770 0.000 0.000
3 -0.565 0.000 0.000 -0.770 0.000 0.000
4 0.565 0.000 0.000 0.770 0.000 0.000
DAV: 8 -0.133293620177E+03 0.15284E-05 -0.29410E-08 4188 0.144E-03 0.119E-04
实际的 DFT 总能量是 OUTCAR 中的能量（例如，energy（sigma->0））减去这个 E_p。根据经验，这是因为 OUTCAR 的能量已经包括了 E_p（我找不到说明这一点的手册，但一些对不同 LAMBDA 的简单测试能够验证这一说法）。我会非常仔细地阅读那本手册，以了解其中的许多微妙之处。
对于 LAMBDA 收敛，我通常会测试一系列值，如 [0，1，10，25，50，75，100，125，150]，并使用合适的值继续进行计算。
根据我的经验，过高的值会导致错误。继续减小范围，例如，减小到 [1,10,25]，然后减小到 [1,3,5,7,10,13,15,17,20]，再减小到 [4,5,6,7,8,9,10,11,12] 等，直到找到最佳值。目的是找到足够大的 LAMBDA，使 E_p 非常小，同时给出所需方向上的最终磁矩。
您可以通过检查 OUTCAR 中 magnetization (x)、magnetization (y) 和 magnetization (z) 的最终结果来验证后者。您需要设置参数 LORBIT=11 ，以便在 OUTCAR 中显示此信息。这些磁化应与 INCAR 中 MAGMOM 和 M_CONSTR 标签定义的数值一致。请注意，如果你不限制大小/幅度，你的幅度可能会有所不同。例如，对于 MAGMOM = 0 0 1 0 1 0 0 1 24*0，我得到：
magnetization (z)
# of ion s p d tot
------------------------------------------
1 0.002 0.002 0.212 0.217
2 0.002 0.002 0.213 0.217
3 0.002 0.002 0.212 0.216
4 -0.000 -0.001 0.001 0.001
5 -0.000 -0.001 0.001 0.001
.....
magnetization (x) 和 magnetization (y) 中的所有离子都为 0。在这种情况下，0.21 的量级是合理的，因为我们预计在这种情况（对于这种特定的材料系统）下，所有 \(\rm Nb\) 原子共享一个半自旋。
这意味着每个 \(\rm Nb\) 原子的大小应该大致为 \(\displaystyle \frac{1}{2} / 3=\frac{1}{6} \approx 0.17\)。
通过改变 RWIGS，可以使这个 0.21 非常接近 0.17。在这种情况下，RWIGS 必须做得更小，这样集成的总面积就会更小，并且不会重叠。但请小心使用这个规律，因为在设置 I_CONSTRAINED_M 后的 OSZICAR 输出页面中还有其他几个相关参数 M_int 和 MW_int。但是，如果使用 I_CONSTRAINED_M=2（也约束大小），则预期最终的力矩大小接近你在 MAGMOM 和 M_CONSTR 中定义的磁矩大小。
最后，您现在拥有了进行约束磁矩计算所需的所有成分。请注意，磁矩收敛过程中可能会有很多问题，互联网有许多潜在的解决方案。这可能是一个棘手的问题，正如 VASP 手册上的一页所示，“祈祷”是一种潜在的解决方法。2
当系统地测试不同的矩方向时，我发现一种非常有用的策略是首先进行自旋非极化计算 ISPIN=1，然后使用其 CHGCAR 和 ICHARG=1 从该电荷密度开始。它似乎有助于收敛、速度和准确性。但是，在使用此 CHGCAR 之前，请确保您的非极化计算使用 OUTCAR 文件检测到预期的对称性。
最后，我想以一种突然的方式分享一个简单的 Python 脚本，我用它来帮助我为 INCAR 文件写出 MAGMOM 和 M_CONSTR 标签。对于大型系统，手动书写可能很麻烦（也很危险！）。所以，希望以下内容有用！
例如 \(\rm Nb_3Cl_8\) 的 \(4\times4\) 的超胞，含有 48 个 \(\rm Nb\) 和 128 个 \(\rm Cl\)。
下面，我给出了前 3 个 \(\rm Nb\) 离子 1 0 0，接下来的 3 个 \(\rm Nb\) 原子 -1 0 0，剩下的 42 个 \(\rm Nb\) 原子 0 1 0，以及所有 \(\rm Cl\) 原子 0 0 0：
def generate_string(a, b, c, n1, n2, n3):
result = []
result.extend(a * n1)
result.extend(b * n2)
result.extend(c * n3)
return "M_CONSTR = " + " ".join(map(str, result)) + " 384*0" # hardcoding 128 * 3 = 384 for Cl
n1 = 3; n2 = 3; n3 = 42 # number of ions
a = [1, 0, 0]; b = [-1, 0, 0]; c = [0, 1, 0] # desired magnetic moments for Nb
print(generate_string(a, b, c, n1, n2, n3))