VASP中如何固定磁矩大小和方向

为了研究材料体系中的相互作用，我们需要研究一些假体系，例如，明明是有自旋轨道耦合（Spin-orbit coupling，SOC）相互作用的Bi2Se3我们会先不加SOC再加SOC，研究自旋轨道耦合对能带反转的作用。在有些体系中我们还会用到调节自旋轨道耦合强度的来研究自旋轨道耦合强度对体系性质的影响。

天玑算老王本王，公众号：科研那些事儿 VASP中只考虑部分元素的自旋轨道相互作用

对于一个给定材料体系的磁矩而言，如果 不考虑外场的影响，只考虑材料最初始/稳定的状态， 可以说一个材料组分和结构确定后就应该只有一种稳定的磁构型（磁矩大小和方向），即磁基态。然而，物理人喜欢研究不同的相之间的相变，因此，计算过程中各种磁构型都可能会计算。

但是对于第一性原理软件VASP而言，只设置ISPIN=2，LSORBIT=T，MAGMOM= 3*Nions*moment。磁矩的大小和方向都有可能被优化到基态。如何实现固定磁矩方向和大小就成了一个问题。

于是，限制磁矩大小和方向的方法被开发出来。见：

https://www.vasp.at/wiki/index.php/Constraining_local_magnetic_moments

VASP中给的案例是Fe的dimer，

POSCAR：

Fe dimer

1.00000000000000000

8.0000000000000000 0.0000000000000000 0.0000000000000000

0.0000000000000000 8.0000000000000000 0.0000000000000000

0.0000000000000000 0.0000000000000000 8.0000000000000000

2

Cartesian

3.00 0.00 0.00

5.00 0.00 0.00

INCAR：

SYSTEM= Fe dimer

ISTART= 0

ISYM= 0

LNonCOLLINEAR = .TRUE.

MAGMOM= 0 0 3 0 0 3

VOSKOWN = 1

LORBIT= 11

! mix slowly when increasing LAMBDA

# AMIX= 0.1

# BMIX= 0.00001

# AMIX_MAG = 0.2

# BMIX_MAG = 0.00001

# I_CONSTRAINED_M = 1

# RWIGS = 1.0

# LAMBDA = 10

# M_ConSTR= 0 0 1 0 0 1

注意，体系设置ISYM=0关了对称性，然后设置两个Fe原子具有001方向的3μB 初始磁矩。

KPOINTS：

k-points

0

Monkhorst Pack

1 1 1

0. 0. 0.

利用MAGMOM = 0 0 3 0 0 3的初始初始铁磁态，在OSZICAR文件中得到如下磁矩:

...

DAV: 20 -0.929676054634E+01 -0.26101E-03 -0.16780E-03 60 0.102E-01 0.537E-02

DAV: 21 -0.929679955346E+01 -0.39007E-04 -0.30319E-04 60 0.590E-02

1 F= -.92967996E+01 E0= -.93047629E+01 d E =0.238900E-01 mag= -0.0006 -0.0003 6.0537

利用MAGMOM = 0 0 3 0 2 2的初始初始铁磁态，在OSZICAR文件中得到如下磁矩:

magnetization (y) magnetization (z)

# of ion s p d tot # of ion s p d tot

---------------------------------------- ----------------------------------------

10.018 -0.001 1.071 1.087 10.045 -0.003 2.587 2.628

20.019 -0.001 1.069 1.087 20.045 -0.003 2.588 2.629

---------------------------------------- ----------------------------------------

tot0.037 -0.003 2.140 2.174 tot0.089 -0.007 5.175 5.257

在系统中加入一个惩罚函数，将体系的局部磁矩驱动到期望的方向，当在输入中修改以下步骤时(注意惩罚函数对总能量的贡献):

磁矩约束的开关(I_CONSTRAINED_M=1)。

设置积分半径来确定局部矩(RWIGS=1.0)。

惩罚函数中的权重(LAMBDA=10)。

磁矩约束的目标方向(M_ConSTR= 0 0 1 0 1 1)。

在OSZICAR中的关键信息为：

E_p = 0.35424E-02 lambda = 0.100E+02

ionMW_int M_int

1 0.000 0.013 1.557 0.000 0.014 2.674

2 0.000 1.092 1.110 0.000 1.880 1.901

DAV: 35 -0.905322335169E+01 0.58398E-04 -0.60872E-04 60 0.734E-02

1 F= -.90532234E+01 E0= -.90355617E+01 d E =-.529849E-01 mag= -0.0005 2.1161 5.1088

E_p是惩罚函数产生的能量。随着LAMBDA增加而减小。

通过逐步增加LAMBDA，可以降低E_p (缓慢地使结果保持稳定):

E_p = 0.22591E-03 lambda = 0.500E+02

ionMW_int M_int

1 0.000 0.002 1.545 0.001 -0.005 2.654

2 0.000 1.086 1.087 0.001 1.871 1.862

DAV: 33 -0.907152551238E+01 0.48186E-04 -0.33125E-04 60 0.163E-01

1 F= -.90715255E+01 E0= -.90541505E+01 d E =-.521251E-01 mag= 0.0042 2.0902 5.0659

这样就可以得到给定磁结构的LSDA总能量。

当收敛性不好时，我们能做什么:

1. 从ISTART=0 (或删除WAVECAR文件)和ICHARG=1计算非自旋极化的电荷密度开始。

2. 通过设置BMIX=0.0001和BMIX_MAG=0.0001来使用线性混合。

3. 慢慢混合，即减少AMIX和AMIX_MAG。

4. 减少MAXMIX，存储在Broyden混合器中的步骤数(默认MAXMIX=45)。

5. 从部分收敛的结果重新启动(大约20步后停止计算，然后从WAVECAR文件重新启动计算)。

6. 使用约束来稳定磁配置。

注意：

当I_CONSTRAINED_M=1时，只固定磁矩方向不固定大小， M_ConSTR= 0 0 1 0 1 1表示方向。

当 I_ConSTRAI NED_M=2时，既固定磁矩方向也固定磁矩大小， M_ConSTR= 0 0 3 0 0 3表示方向和大小。