自旋极化和旋轨耦合相关问题

asl1994

初学VASP，目前想重现一个文献中对PuO2的DOS和能带结构的计算，但是有些基本的问题还没有弄清楚。1）在PuO2里Pu是+4价而O是-2价，而Pu原子价电子为5f(6)4s(2)，Pu4+就是5f(4)吧。可是为什么我设置LORBIT=11在OUTCAR结果看到如下：
total charge

# of ion     s       p       d       f       tot
------------------------------------------------
  1        2.180   5.982   1.142   4.852  14.156
  2        2.180   5.982   1.142   4.852  14.156
  3        2.180   5.982   1.142   4.851  14.155
  4        2.180   5.982   1.142   4.852  14.155
  5        1.560   3.457   0.000   0.000   5.017
  6        1.560   3.457   0.000   0.000   5.017
  7        1.560   3.457   0.000   0.000   5.017
  8        1.560   3.457   0.000   0.000   5.017
  9        1.560   3.457   0.000   0.000   5.017
10        1.560   3.457   0.000   0.000   5.017
11        1.560   3.457   0.000   0.000   5.017
12        1.560   3.457   0.000   0.000   5.017
------------------------------------------------
tot       21.200  51.585   4.568  19.406  96.759


这里Pu的价电子要怎么看呢？似乎和推测的不符啊。
赝势开头是
  ##PAW Pu 17Apr2000
   16.0000000000000
parameters from PSCTR are:
   VRHFIN =Pu:  [Xe,5d,4f]
   LEXCH  = CA
   EATOM  =  1953.6839 eV,  143.5917 Ry

这样的赝势有问题吗？

2）按照图1的说法，Pu4+的4个5f电子都成对，PuO2作为抗磁（diamagnetic）体系来计算。而作者又说在计算里都考虑了自旋-轨道相互作用。我想知道既然没有单电子，计算需要开启自旋吗？我试着开了自旋（ISPIN=2），不设置初始磁矩，算出来的总磁矩有16μB（如下）
magnetization (x)

# of ion     s       p       d       f       tot
------------------------------------------------
  1       -0.014   0.000  -0.063  -4.323  -4.399
  2       -0.014   0.000  -0.063  -4.323  -4.399
  3       -0.014   0.000  -0.063  -4.322  -4.399
  4       -0.014   0.000  -0.063  -4.322  -4.399
  5        0.007   0.158   0.000   0.000   0.165
  6        0.007   0.158   0.000   0.000   0.165
  7        0.007   0.158   0.000   0.000   0.165
  8        0.007   0.158   0.000   0.000   0.165
  9        0.007   0.158   0.000   0.000   0.165
10        0.007   0.158   0.000   0.000   0.165
11        0.007   0.158   0.000   0.000   0.165
12        0.007   0.158   0.000   0.000   0.165
------------------------------------------------
tot       -0.001   1.266  -0.252 -17.290 -16.276

这是意味着每个Pu上有4个未成对的f电子吗？如果再计算旋轨耦合结果会像图1一样吗？

静态自洽的INCAR如下：
SYSTEM = PuO2

ISTART = 0
ICHARG = 2
ISPIN = 2
VOSKOWN = 1
MAGMOM = 8*0 4*0
LORBIT = 11
LREAL = .F.
NCORE = 4
LSORBIT = .F.

# Accuracy controls:
PREC = normal
ENCUT = 500

# Electronic loop:
ALGO = Normal    # Normal (Davidson)
EDIFF = 1E-07
NELM = 100     # Maximum number of electronic SC steps (default:60)

# DOS related values:
ISMEAR = -5        # tetrahedron method with Blochl corrections (use a G-centered k-mesh)
SIGMA = 0.200000   # eV

# Ironic relaxation:
IBRION = -1
ISIF = 2       # Cell fixed
NSW = 0 # Maximum number of ionic steps
POTIM = 0.500000 (scaling constant)
EDIFFG= -0.01

# DFT +U
LDAU = .T.
LDAUTYPE = 2
LDAUL = 3 -1
LDAUU = 6.35 0
LDAUJ = 0 0

# Write flags:
LCHARG = .True.
LWAVE = .True.
LELF = .F.

LVTOT = .F.
LVHAR = .F.


然后尝试加SOC算DOS：

ISTART = 1
ICHARG = 11
ISPIN = 2
VOSKOWN = 1
LORBIT = 10
LREAL = .F.
NCORE = 4
LSORBIT = .True.
LMAXMIX = 6
ISYM = 0
SAXIS = 0 0 1

# Accuracy controls:
PREC = normal
ENCUT = 500

# Electronic loop:
ALGO = Normal    # Normal (Davidson)
EDIFF = 1E-07
NELM = 100     # Maximum number of electronic SC steps (default:60)

# DOS related values:
ISMEAR = -5        # tetrahedron method with Blochl corrections (use a G-centered k-mesh)
SIGMA = 0.200000   # eV
NBANDS = 140
EMIN = -50
EMAX = 15
NEDOS = 2000

# Ironic relaxation:
IBRION = -1
ISIF = 2       # Cell fixed
NSW = 0 # Maximum number of ionic steps
POTIM = 0.500000 (scaling constant)
EDIFFG= -0.01

# DFT +U
LDAU = .T.
LDAUTYPE = 2
LDAUL = 3 -1
LDAUU = 6.35 0
LDAUJ = 0 0

但是结果马上出错了。
POSCAR found type information on POSCAR  Pu O
POSCAR found :  2 types and      12 ions
scaLAPACK will be used
LDA part: xc-table for Ceperly-Alder, Vosko type interpolation para-ferro
POSCAR found type information on POSCAR  Pu O
POSCAR found :  2 types and      12 ions
found WAVECAR, reading the header
  number of k-points has changed, file:    35 present:   729
  trying to continue reading WAVECAR, but it might fail
POSCAR, INCAR and KPOINTS ok, starting setup
WARNING: small aliasing (wrap around) errors must be expected
FFT: planning ...
reading WAVECAR
reading wavefunctions of collinear run, up
reading wavefunctions of collinear run, down
the WAVECAR file was read successfully
charge-density read from file: PuO2
magnetization density read from file 1
LAPACK: Routine ZPOTRF failed!           1          36           1
LAPACK: Routine ZPOTRF failed!
LAPACK: Routine ZPOTRF failed!           1          36           1
LAPACK: Routine ZPOTRF failed!
LAPACK: Routine ZPOTRF failed!           1          36           1
LAPACK: Routine ZPOTRF failed!
LAPACK: Routine ZPOTRF failed!           1          36           1
LAPACK: Routine ZPOTRF failed!
LAPACK: Routine ZPOTRF failed!           1          36           1
LAPACK: Routine ZPOTRF failed!
LAPACK: Routine ZPOTRF failed!           1          36           1
LAPACK: Routine ZPOTRF failed!
LAPACK: Routine ZPOTRF failed!           1          36           1
LAPACK: Routine ZPOTRF failed!
LAPACK: Routine ZPOTRF failed!           1          36           1
LAPACK: Routine ZPOTRF failed!
LAPACK: Routine ZPOTRF failed!           1          36           1
LAPACK: Routine ZPOTRF failed!

请问是哪里出了问题呢？

丁越

卡开发发 发表于 2022-1-28 18:22
没有没有，说不上指导，不懂的很多，欢迎讨论。这方面的原理你想进一步了解可以看看Richard M Martin的《 ...

老师太谦虚了哈哈哈，非常感谢

卡开发发

丁越 发表于 2022-1-28 18:10
感谢卡大人指导

没有没有，说不上指导，不懂的很多，欢迎讨论。这方面的原理你想进一步了解可以看看Richard M Martin的《电子结构》，第11章有模型势/赝势/PAW专门的内容，但是提及的不够全面，书后的文献可以补足这些内容。我自己之前因为用过一段时间SIESTA的ATOM生成Troullier-Martins模守恒赝势，找过一些资料研究了一下。

丁越

卡开发发 发表于 2022-1-28 16:12
这个需要测试，如果对你的体系是否包含半芯态对前线电子的电子结构影响很小，那么没有必要去考虑半芯态； ...

感谢卡大人指导

卡开发发

丁越 发表于 2022-1-28 14:18
请教一下卡开发发老师，包含半芯态电子的赝势和不包含半芯态电子赝势在计算中该怎么选择呢（假设是QE中的 ...

这个需要测试，如果对你的体系是否包含半芯态对前线电子的电子结构影响很小，那么没有必要去考虑半芯态；反之，我建议还是和全电子结果小心比对再考虑使用什么。

半芯态就是教材里面提到的“小核势”对应，也就是说增设了内层电子，一般相应而言芯半径也会小一些，动能截断也会大一些以适应更复杂的化学环境。“大核势”虽然牵涉更少的电子，但对于不同化学环境的可迁移性没那么好，也就是说特定的化合物（或者说特定的化合价与杂化方式）中，这些赝势不能复现与全电子一致的结果。通常对于形状一致赝势，生成赝势后其实会对不同的原子价层排布来验证这种可迁移性，但是验证的电子排布方式是有限的，而且也只是针对孤立的原子情形。

QE的赝势应该还是经过delta test检验过了的，但我也不能保证即便delta test还行的赝势到了某个特定体系表现一定能不错。另外要小心的就是，我个人测试下来，QE某些赝势还是要对空带数目做些调整。总体来说小心测试一定是对的，尤其是那些开源的程序。

丁越

卡开发发 发表于 2017-9-21 12:49
1、PAW的组态没问题。
（1）一般过渡元素为了同时保证PAW的精度和可移植性，会把半芯态也作为active elect ...

请教一下卡开发发老师，包含半芯态电子的赝势和不包含半芯态电子赝势在计算中该怎么选择呢（假设是QE中的Fe元素的PAW赝势，Fe.pbe-n-kjpaw_psl.1.0.0.UPF和Fe.pbe-spn-kjpaw_psl.1.0.0.UPF）？即哪些场景下必须得用包含半芯态的赝势，而哪些场景下计算时不包含半芯态的赝势就足够了？

卡开发发

asl1994 发表于 2017-9-21 16:18
谢谢老师！按照您说的我先用ISYM=0做了自洽计算，初始磁矩设为零，算完发现总磁矩又变成2了，OUTCAR结果 ...

1、看起来好像不是很正常，难道计算完之后四个Pu的位置会出现不等价？实在不行你可以ISPIN=1做，然后SOC的时候再改过来。
2、EIGENVAL显示的这个问题我不太清楚是怎么发生的，实在不行也只好到vasp的社区去问一下。

asl1994

卡开发发 发表于 2017-9-21 12:49
1、PAW的组态没问题。
（1）一般过渡元素为了同时保证PAW的精度和可移植性，会把半芯态也作为active elect ...

谢谢老师！按照您说的我先用ISYM=0做了自洽计算，初始磁矩设为零，算完发现总磁矩又变成2了，OUTCAR结果对应的是 magnetization (x)

# of ion     s       p       d       f       tot
------------------------------------------------
  1        0.008   0.005   0.025   2.012   2.050
  2       -0.007  -0.006  -0.024  -2.042  -2.079
  3        0.001  -0.000   0.003  -0.001   0.003
  4        0.007   0.005   0.030   2.062   2.104
  5       -0.001  -0.011   0.000   0.000  -0.013
  6       -0.001  -0.011   0.000   0.000  -0.013
  7       -0.001  -0.012   0.000   0.000  -0.013
  8       -0.001  -0.012   0.000   0.000  -0.013
  9       -0.000  -0.003   0.000   0.000  -0.003
10       -0.000  -0.003   0.000   0.000  -0.003
11       -0.000   0.001   0.000   0.000   0.000
12       -0.000   0.001   0.000   0.000   0.000
------------------------------------------------
tot        0.004  -0.048   0.034   2.031   2.020

这个结果正常吗？

接着我把自洽计算的IBZKPT作为非自洽算SOC的KPOINTS，计算正常结束了，可是总磁矩还是2。要怎么做才能得到文献中所有电子都成对的情况呢？
另外而且SOC算完EIGENVAL里只有三列，最后一列占据数是1，这是意味着已经没有alpha和beta自旋的区别了吗？
还有一点很奇怪，我最开始用ISPIN=1算，自洽完的EIGENVAL是三列但最后一列也是1？不是应该是2吗？

ISPIN=1自洽完的EIGENVAL前几行：
   12   12    1    1
  0.1310881E+02  0.5398190E-09  0.5398190E-09  0.5398190E-09  0.5000000E-15
  1.000000000000000E-004
  CAR
PuO2
    112     35     70

  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.1371742E-02
    1      -37.075216   1.000000
    2      -37.002301   1.000000
    3      -37.002301   1.000000
    4      -37.002301   1.000000
    5      -12.801582   1.000000
    6      -12.801582   1.000000
    7      -12.801582   1.000000
    8      -10.204648   1.000000
    9      -10.204648   1.000000
   10      -10.204648   1.000000

（为什么占据数是1？）

ISPIN=2自洽完的EIGENVAL前几行：   12   12    1    2
  0.1310881E+02  0.5398190E-09  0.5398190E-09  0.5398190E-09  0.5000000E-15
  1.000000000000000E-004
  CAR
PuO2
    112    365     70

  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.1371742E-02
    1      -38.507498    -38.268304   1.000000   1.000000
    2      -38.024082    -37.748500   1.000000   1.000000
    3      -37.754633    -37.558535   1.000000   1.000000
    4      -37.267291    -36.950175   1.000000   1.000000
    5      -13.091142    -12.990377   1.000000   1.000000
    6      -13.005922    -12.856144   1.000000   1.000000
    7      -12.868493    -12.798521   1.000000   1.000000
    8      -11.268173    -10.954608   1.000000   1.000000
    9      -11.058985    -10.898223   1.000000   1.000000
   10      -11.031836    -10.811231   1.000000   1.000000


ISPIN=2算SOC的EIGENVAL前几行：
   12   12    1    1
  0.1310881E+02  0.5398190E-09  0.5398190E-09  0.5398190E-09  0.5000000E-15
  1.000000000000000E-004
  CAR
PuO2
    112    365    140

  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.1371742E-02
    1      -38.515293   1.000000
    2      -38.279267   1.000000
    3      -38.048236   1.000000
    4      -37.778721   1.000000
    5      -37.758864   1.000000
    6      -37.557159   1.000000
    7      -37.274320   1.000000
    8      -36.970548   1.000000
    9      -17.898483   1.000000
   10      -17.679977   1.000000

卡开发发

1、PAW的组态没问题。
（1）一般过渡元素为了同时保证PAW的精度和可移植性，会把半芯态也作为active electrons（这个概念就是实际参与计算的电子），这个active electrons并不等于我们一般定义的价电子数目就是这个原因。有时候生成赝势/PAW的电子组态也不一定是按照原子基态产生的，比如Fe的基态是3d6 4s2，但是实际PAW/赝势是按照3d7 4s1的组态产生。
（2）做布居分析也不会恰好和价态一样，不同布居分析的处理方法结果数值上也基本不同。

2、（1）你如果不设置磁矩，对应ISPIN=2默认就是铁磁性的

1. Default: MAGMOM        = NIONS * 1.0

复制代码

你应该计算铁磁性和无磁性的情况，然后比较SCF能量。
（2）通过非自洽计算SOC之前，自洽计算的时候尝试ISYM=0试试看看。