7.1 磁性计算1

发表于 更新于 分类于

vasp磁性参数

设置ISPIN和MAGMOM参数，详见vasp官网。关闭对称性：ISYM=0.一般有的采用先做非磁结构优化，再打开自旋和磁性设置做磁基态计算。确定磁性体系可以直接做磁性计算。磁性计算一般属于过渡金属，要进行DFT+U计算。对于考虑LDAU的情况一般推荐增加参数 LMAXMIX。

INCAR中的参数设置：

1
2
3
4
5
6
7
LORBIT = 11
LDAU = .TRUE. 打开LDAU
LDAUTYPE =2 LDAU类型LDAUTYPE=2： 只考虑U-J的值；LDAUTYPE=1： 同时考虑U和J，J基本为常数1 eV
LDAUL = 2 -1 库伦排斥的轨道
LDAUU = 2.8 0 库伦排斥的大小
LDAUJ = 0 0 stoner 交换参数大小
LMAXMIX = 4 mix轨道

如果进行非线性计算，LNONCOLLINEAR = .TRUE. MAGMOM = x y z . . . … 可以读前一步的CHGCAR，进行非线性自洽计算（不用结构优化），也可以直接进行非线性计算。为了避免以上报错，增
加mix参数：

1
2
3
4
5
LMAXMIX = 4
AMIX = 0.2
BMIX = 0.00001
AMIX_MAG = 0.8
BMIX_MAG = 0.00001

提交任务要使用vasp_ncl版本的程序

自旋轨道耦合
默认 SAXIS = 0 0 1
加入磁矩参数
MAGMOM = x y z … …

一般接口wannier90还需要
ISYM=-1
NPAR =1 或者将这个开关删除
##

铁磁半导体相对好收敛，遇到不收敛问题，可尝试降低 AMIX，增加 BMIX
尝试更换不同的ISMEAR，比如加电场时，会遇到矩阵不厄密。
call to ZHEGV failed 问题可能是

  1. 结构不合理
  2. 优化时优化到一个不合理的结构： 尝试改变弛豫算法IBRION ，减小步长POTIM = 0.1或者更小
  3. 对角化算法不稳定 ： ALGO = Normal/Fast
  4. 尝试关掉对称性

磁基态

得出体系的具体磁基态：比较三个计算无磁，铁磁，反铁磁的能量，能量更低的为体系的磁基态。构建不同磁基态有可能需要扩胞。

得出体系原子的具体磁矩：用vi编辑器打开OUTCAR或其他方式打开，从最后往前翻动，第一个magnetization (x)开头的表格就是原子的磁矩大小。需要打开LORBIT=11

磁矩分布画图：使用chgsplit.sh（可以向组内师兄师姐寻找）。磁性测试输出CHGCAR后，运行命令./chgsplit.sh CHGCAR得到cf1和cf2文件，把cf2文件重命名为CHGCAR-mag，与CONTCAR一起拖入VESTA中画图即可。

磁耦合

对于磁性体系的第一性原理计算结果，能量部分可以分
为两个部分E = E_0 （非磁项） + E_s （磁性相关能量）
对比磁性耦合产生的能量，不能够用磁性态的能量跟非
磁的去对比
E_s 包含交换耦合以及单离子各向异性能（SOC）
一般来说，磁性相关的能量可以用下面的公式来表示

  1. 磁性态的构造
    为了求磁耦合强度，我们需要构造出不同磁构型哈密顿，通过求解这些哈密顿方程
    组来得到磁耦合参数。例如对于一个体系，我们需要求解最近邻J1和次
    近邻J2，那么体系能量中有三个未知量H0, J1, J2。这样我们需要构造三个线性独立的方程来求解。
    磁耦合参数

四态法

接下来，我们以一个例子来看四态法
是如何用来求磁耦合的。

  1. 构建超胞
    这里我们以二维六角晶格为例，原胞
    中有两个磁性原子，我们构建一个
    881的超胞，目标是求下图中A、B
    之间的磁耦合作用。
    构建四个磁性态 ，分别
    为
    其中第一项是DFT中的非磁项，联立这四个方程组就
    可以求解
    注意这里的超胞大小，只要满足这个8*8的超胞里的A、
    B与周期外的A、B格点之间没有作用或者作用可忽略就
    可以。

单离子各向异性，对于沿Z轴具有三重四重六重旋转对称性的体系，只需要计算一个Azz − Axx 。

居里温度

对于铁磁有序材料，当温度达到某一临界值之后，在没有外磁场的情况下，材料会变成
顺磁态，这个温度点便是居里温度（Curie temperature: Tc）
铁磁是一种铁性材料，所谓铁磁是指在外场作用下，极化随着外场的翻转会有一个滞后，
自旋极化随着外场的变化曲线形成磁滞回线。对于反铁磁有序材料，当温度达到某一临界值之后，在没有外磁场的情况下，材料同样会变成
顺磁态，这个温度点便是奈尔温度（Neel temperature: TN）

对于反铁磁有序材料，一般可以使用测试磁化率来判断其相变温度，磁化率满足居里外斯定律

  1. 平均场近似
  2. 自旋波方法详细推导见文件夹mean-field-theory
  3. Monte Carlo方法

居里温度的Monte Carlo 模拟

  1. 撒点，给初始磁矩
  2. 计算哈密顿
    3 尝试随机反转一个格点的磁矩，计算能量
    4判断，如果能量变低，则翻转，如果能量变高，那么
    如果 则翻转，否则不翻转
    5 先进行充足的MC步骤至平衡，然后在进行充足的步骤采样，取平均磁矩。
    每一个温度点做一次1-5步骤，即每一个温度点有一个平均磁矩，磁矩下降最快的点（斜率绝对值最大）即为
    居里温度。
    反铁磁需要计算奈尔温度，那么平均磁矩会一直是0，有两种办法，算其中一个格点的平均磁矩，或者看磁比
    热

mcsolver计算XY模型的一个例子。

在学习的过程中，我发现参数设置非常重要，也比较难，但网上并没有给出一个详细的教程说明，结合网上已有说明和朋友讨论，以下给出我的一点见解，希望大家相互讨论学习，有错误欢迎批评指正。

1.lattice晶格常数

关于晶格常数的设置，这步比较简单。

晶格常数需要采用归一化后的晶格常数，a1、a2、a3分别归一化，例如a1=（3，4，0）归一化后就是（0.6，0.8，0）。

sc（supercell）表示使用的超胞规模，一般情况下尽可能的大，例如16161或32321（在二维情况下）。

2.Orbital list

这一步是对原胞中不等价磁性原子的设置。

ID是对设置的原子的编号，在可视化界面中，这一项不用手动设置，程序会自动根据我们的设置进行从0开始的编号。

type表示原胞中磁性原子的种类，根据具体情况进行设置，如果只有一种磁性原子，都设置成0就可以。

init spin表示自旋磁矩的设置，比如你用VASP计算出来得到的磁矩是3，那么这一项需要设置为1.5（需要除以一个2，换算成玻尔磁子）。

pos对磁性原子的位置坐标进行设置，这里最好是采用分数坐标（如果是笛卡尔坐标建议转换成分数坐标进行计算）。（更正：应该是相对于a1、a2、a3的分数坐标，也就是相对于斜的平行四边形的而言的，以其两条邻边为晶格矢量，例如CrI3是0.333 0.666667 0，如果是直角坐标系下的分数坐标的话第一个值应该是负的才对）

Ani表示单离子在xyz方向上的各向异性系数或单离子磁晶各向异性能(在Ising模型中是无用的，在XY模型中只使用前两个）。对自己的结构做一个易磁化轴位置的判断。一般为Z方向就在DZ处添加。同时，注意各向异性的单位是开尔文。 （1meV=11.604609K）

至此，我们已经完成了第二步Orbital list的设置。

3.Bond list

设置第二步里面的磁性原子之间的交换耦合系数，在你的计算中，一共考虑到了多少个相互作用就要添加多少个。比如，对于CrI3来说，考虑到了最近邻、次近邻和第三近邻相互作用，分别有3个、6个和3个，那么你一共就需要添加12条信息。

对于每一项来说，设置方法如下：

ID是相互作用的编号，如同第二步一样，系统会从0开始自动编号。

J是磁性原子之间的交换耦合系数，一个J有九个矩阵元素，分别包括Jxx、Jyy、Jzz、Jxy、Jxz、Jyz、Jyx、Jzx、Jzy。每个元素描述了自旋的两个分量之间的耦合。对于Ising模型，由于只考虑一维自选变数，即只存在自旋向上或自旋向下，因此只使用第一个元素Jxx。对于XY模型也一样，自旋的朝向由up和down解放到了XY平面的任意朝向，因此只使用Jxx、Jyy、Jxy、Jyz。而对于heisenberg模型来说，三个自旋方向都存在相互作用，因此Jxx、Jyy、Jzz、Jxy、Jxz、Jyz、Jyx、Jzx、Jzy九个矩阵元素都需要设置。（对于模型的设置在后面的Model处进行选择。）对于前Ising和XY模型来说，即便你设置了九个矩阵元素，有效输入的也只是对应的起作用的元素，因此选好模型很关键。

s，t和over lat.表示我们考虑哪两个磁性原子之间的相互作用，s和t要选择我们在第二步中设置的原胞中磁性原子的ID，over lat.的三个元素指代一个晶格矢量，这个晶格矢量表示：在晶格中，你考虑的原子（在该中心原胞平移该晶格矢量后的原胞的t原子），和中心s原子之间的相互作用。（这里比较难，为了便于理解，建议读者找个例子琢磨琢磨或者自己多动手尝试一下）

注意所有的能量单位都是开尔文（1meV=11.604609K）

至此，我们已经完成了第三步Bond list的设置。

现在，我们可以在Structure viewer里面看到超胞里磁性原子的分布和相互作用示意图。大家可以根据右上角的图查看自己添加的是否正确。一定要多加尝试。

4.Other settings其他设置

关于一些其他参数的设置。

T start&end表示起始和结束温度，也就是温度的取值范围，total points表示温度插值次数（用于温度扫描），也就是总取点数，一般来说点越多图形越精确。示例给出的0.9,1.2和8只是为了快速计算得到结果，真实输入的时候要视实际情况而定。

同理，H start&end表示外加磁场的取值范围，total points表示采样次数（用于磁场扫描）。

nthermal是使系统进入平衡状态的总步骤，nsweep是测量所涉及的总步骤，tau表示每一步的MC更新。

xAxis表示右边Result viewer的x轴上的物理量，可以是T(表示M-T曲线)或H(表示迟滞回线)。

Model可以选择Ising模型、XY模型或Heisenberg模型。

Algorithm选择相应的算法（支持Wolff，Metropolis，Swedsen-Wang）。

Measure corr. si&sj设置spin＿i和spin＿j以及它们之间的晶格矢量，用于相关拼接。（如果spin＿i=spin＿j并且overLat=0 0 0，那么你将得到spin＿i的磁化率）。

nFrame是输出自旋构型的数目，用于说明在平衡或非平衡状态下的自旋构型。

core设置并行计算的核心资源。

至此，所有参数都设置完毕。

5.最后

save可以选择保存当前参数到文件，以便下一次找到。

单击StartMC启动计算。

等待右面板中的关系图更新。之后，你可以在mcsolver的根目录下找到一个result.txt文件，其中有许多有用的信息，包括平均自旋(在步骤5中定义的spin＿i和j上)，spin＿i和j之间的相关性，内能，比热容量和Binder累积量U4等。如果你处理有多个核心的模拟，那么结果可能不会根据温度排序，但是，每一行的对应关系都是可以的。

原文链接：https://blog.csdn.net/qq_46679797/article/details/134455309

MAE

磁晶各向异性 magnetic anisotropic energy (MAE) 是指自旋方向在不同方向上的能量差。MAE跟体系对称性有关。
对于简单体二维系，易磁化轴可能在面内或者垂直于面，我们可以选择计算不同方向的某一磁性态能量对比来求MAE
比如CrI3， 我们可以计算沿着x方向的FM能量Ex，然后计算沿着z方向的FM能量Ez。
MAE = (Ex -Ez)/2
对于对称性较低的体系，需要每个面去计算MAE。在vasp中，三维体系研究MAE，可以在高对称轴上算能量差

三步：结构优化，共线计算磁性并输出CHGCAR,读取CHGCAR改为非线性计算不同方向的能量。

晶体磁各向异性能量由按不同方向旋转所有自旋决定。首先，必须在基态下进行精确（PREC=精确，LREAL=.False.）共线计算（使用vasp_std版本）。接下来，考虑自旋轨道耦合（LSORBIT=.True.，使用vasp_ncl版本）对几种自旋取向进行非自洽计算（ICHARG=11）。在大多数情况下，能量变化非常小（有时约为微电子伏）。与共线相比，必须将NBANDS加倍。

为了修改晶体中自旋的取向，我们考虑这里描述的第二种方法。对于MAGMOM标签，根据z方向编写总局部磁矩（必要时，x和y方向均为0）。自旋取向（𝑢，𝑣，𝑤）由笛卡尔坐标系中的SAXIS标签定义。通过使自旋朝不同方向取向计算晶体磁各向异性能量和以下方程式

𝐸MAE=𝐸（𝑢，𝑣，𝑤）−𝐸0

其中 𝐸0 是最稳定自旋方向的能量。更多细节请查看SAXIS和LSORBIT页面。

练习：通过将自旋沿以下路径定向，确定NiO的非自洽方式下的晶体磁各向异性能量：(2,2,2) –> (2,2,1) –> (2,2,0) –> … –> (2,2,-6)。与自洽方法进行比较。

INCAR:

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
NiO MAE
  SYSTEM    = "NiO"
    
Electronic minimization
  PREC = Accurate
  ENCUT         = 450
  EDIFF         = 1E-7
  LORBIT        = 11
  LREAL         = .False.
  ISYM          = -1
  NELMIN        = 6
  #  ICHARG = 11
  #  LCHARG = .FALSE.
  #  LWAVE = .FALSE.
  #  NBANDS = 52
  #  GGA_COMPAT = .FALSE.
    
DOS
  ISMEAR    = -5
    
Magnetism
  ISPIN     = 2
  MAGMOM = 2.0 -2.0 2*0.0
  # MAGMOM    = 0 0 2 0 0 -2 6*0 # Including Spin-orbit
  # LSORBIT       = .True.
  # SAXIS = 1 0 0 # Quantization axis used to rotate all spins in a direction defined in the (O,x,y,z) Cartesian frame
    
Orbital mom.
  LORBMOM = T
    
Mixer
  AMIX      = 0.2
  BMIX      = 0.00001
  AMIX_MAG  = 0.8
  BMIX_MAG  = 0.00001
    
GGA+U
  LDAU      = .TRUE.
  LDAUTYPE  = 2
  LDAUL     = 2 -1
  LDAUU     = 5.00 0.00
  LDAUJ     = 0.00 0.00
  LDAUPRINT = 1
  LMAXMIX   = 4

vaspkit

先执行结构优化计算（略），之后使用以下INCAR参数并选择991 K点进行自旋极化静态计算，运行vasp_std得到CHGCAR文件。

MAE的INCAR：

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
#General
PREC = ACCURATE
ISTART= 0
ICHARG= 11
ENCUT = 520
EDIFF = 1E-8
EDIFFG = -0.001
LREAL = .F.
NPAR = 4
NSW= 0
IBRION = -1
ISIF = 2
ISMEAR = -5
SIGMA = 0.05
LCHARG = .F.
LWAVE = .F.
POTIM = 0.5
NEDOS = 2001
NBANDS = 144        # 修改为上一步静态计算NBANDS数值的2倍

#Magnetic properties

MAGMOM = 18*0 0 0 4 0 0 4
NELMIN = 6
LORBIT = 11
ISYM = 0
LSORBIT = .True.
LMAXMIX = 4  ! for d-elements increase LMAXMIX to 4, f-elements: LMAXMIX = 6

和VPKIT.in

1
2
3
1                              ! 1为预处理, 2为后处理
0 360 12                       ! Phi, Spherical coordinate system
0 180 12                       ! 0<= theta <=180, 90 degree means in x-y plane

VPKIT.in说明：
第一行：在第二步计算前使用1，第二步计算完后设置为2；

第二行：三维球坐标系下矢量在xy面投影与x轴正方向的夹角方位角φ，由起始角度0°到终止角度360°(0°)，最后是分割角度的个数；

第三行：三维球坐标系下矢量与z轴的夹角θ，由起始角度0°到终止角度180°，最后是分割角度的个数(注意包含180°在内需要＋1，所以是7)；

其中和分别为方位角和极角，其值变化范围分别是从0°到180°和从0°到360°，这里方位角和极角各取12个离散点（共计144个计算作业）。如果计算资源允许，可以选取更小间隔。调用vaspkit-621命令产生如下所示的一系列计算作业。理论上要生成144个计算作业，考虑到一部分相等无需计算，实际上共生成了122个计算作业。预处理时VASPKIT会自动在每个计算作业文件夹里的INCAR后根据公式(1)追加SAXIS =

接下来可使用以下脚本模板批量提交vasp_ncl作业。

1
2
3
4
5
6
7
8
9
10
#!/bin/bash
# just for reference
path=`pwd`
for c in *.*
do
   cd $path/$c
   echo `pwd`
   qsub vasp_job_soc.sh   ! vasp_job_soc.sh为vasp_ncl作业脚本
  cd ${path}
done

待vasp_ncl计算完成后，把INPUT.in文件中的第一行修改2，然后再次运行vaspkit-621得到MAE.dat文件。我们可利用vaspkit/examples/MAE文件夹中的MATLAB脚本mae_3D_plot_matlab.m进行可视化。

  • 如果想得到xy二维平面内的MAE，INPUT.in可设置为:
    1
    2
    3
    1                              ! 1为预处理, 2为后处理
    0 360 12                       ! Phi, Spherical coordinate system
    90 90 1                        ! xy平面内theta始终等于90°
  • 如果想得到xz二维平面内的MAE，INPUT.in可设置为:
    1
    2
    3
    1                              ! 1为预处理, 2为后处理
    0 0 1                      ! Phi, Spherical coordinate system
    0 180 12                       !
    tips:

根据提示我们需要先执行一步sta，然后读取WAVECAR和CHGCAR进行Nocollinear的计算，当读取WAVECAR进行测试的时候，计算是一定会报错的，因为俩步计算电子的自由度是不同。这时候VASP一定会提示：ERROR: while reading WACECAR,plane wave coefficidents changed

那么继续根据官网的操作手册

  • 我们第一步只计算无磁性基态的CHGCAR,第二步去读取CHGCAR进行计算。
    这里只比对了x方向和z方向的结果：发现计算的结果是一模一样的，无法区分。

  • 接下来不读取之前的CHGCAR进行计算，还是依照之前的设置.这时候发现计算的结果是有区分的。方向不同。

  • 还有一种计算方式是读取写入磁基态的CHGCAR，再进行自旋轨道耦合的计算。这样计算的结果是，无论选择量子化轴在哪个方向，计算的结果都是在Z轴方向，但是能量有所差距。但是根据以上俩种结果去推断Cr原子的磁晶各向异性能，发现这种计算的结果跟文献的600-800eV的结果更为接近。 所以姑且认为进行计算是合理的结果。

电子结构分析

磁性一般来自于d、f轨道的半占据电子，下面以d轨道为例介绍电子结构的分析

  1. 磁矩约等于非成对的电子数，但是电子的占据方式并不仅仅由洪特规则来决定。
    以一个Mn2+为例，最外层5个电子，按照洪特规则Mn磁矩应该是5，但是在晶体场的作用下，它有可能呈现高自旋、
    中自旋、低自旋态。

如何查找结构晶体场：

  1. 确定点群 运行vaspkit
  2. 查表
    http://gernot-katzers-spice- pages.com/character_tables/index.html
  3. 判断能及大小