Li金属表面放带电有机分子体系的弛豫计算输入文件设置分享

Eva

本人近期在学习用VASP计算，现在需要将超氧自由基和环戊酮放在金属Li的（110）表面上。关于四个输入文件的准备过程如下。
1. POSCAR
  通过MS建模后，利用VESTA生成POSCAR文件。需要注意的是，（1）由于有机分子比较长，真空层需要设置大一些（我设置了20/AA）；
（2）由于此体系带一个负电，需要设置NELECT，NELECT是价电子数综合，可通过POTCAR确定每种原子的价电子数。带电这意味着此体系在真空方向上会因为截断周期性而产生虚假的镜像相互作用（本应通过金属表面的感应电荷而终止的偶极子电场，穿过真空遇到周期性产生的另一个偶极子，他们会产生静电作用），需要偶极修正。偶极校正的三个关键词：
LDIPOL =  .TRUE.  （开启偶极校正，supercell需要为cubic）
IDIPOL =  3           (The direction for dipole correction，由于只有Z方向有真空层，就在Z方向上进行修正即可)
DIPOL  = 0.5 0.5 0.5 (The position for dipole correction, fraction coordination，偶极修正的原理是打破平动对称性，当电子密度在Rcenter一定距离趋于0时此方法才能得出cinsistent value，否则偶极校正值会对center的位置非常敏感，此时更建议增加晶胞尺寸或者扩大真空层。center位置建议选在晶胞内原子的质心位置。如果不手动设置，程序会先确定某平面内平均电荷密度降至最小值的位置，并通过将该最小值所在平面的垂直晶格矢量的一半加到此平面上，以此计算出电荷分布的中心位置，从计算方法可知，偶极校正对正交晶胞更可靠。
那么问题来了，supercell需要为cubic（POSCAR中的晶格常数a=b=c），但构建真空层本身不是cubic怎么办呢，通过扩胞计算量会不会太大呢？目前的想法是偶极校正并不影响弛豫计算的结构，先舍弃偶极校正试试。
（3）由于此体系带一个单电子，需要设置MAGMOM，需要注意这里的设置顺序和POSCAR的原子排序以及坐标排序要相互对应。

2. INCAR
（1）弛豫计算会先进行电荷密度初猜，进入SCF循环时先利用电荷密度构建哈密顿矩阵，再对角化得到波函数，接着利用波函数计算新电荷密度，混合新旧电荷密度后再构建哈密顿矩阵，直至达到力收敛标准。ICHARG=0 Calculate the charge density from initial wave functions；ICHARG=1 Read the charge density from CHGCAR file；ICHARG=2 Take superposition of atomic charge densities。 此体系的电荷和自旋多重度需要单独定义，个人目前认为直接算出电荷比较好，设置为ICHARG=2.
（2）
INCAR内容如下：

Global Parameters
ISTART =  0            (Read existing wavefunction, if there)
ISPIN  =  2            (Spin polarised DFT)
NELECT =   167
NUPDOWN = 1
MAGMOM =  40*0.01 1.0 -1.0 3*0.01 3*0.01 4*0.01
ICHARG =  2            (Non-self-consistent: GGA/LDA band structures)
LREAL  = Auto       (Projection ope，rators: automatic)
ENCUT  =  520          (Cut-off energy for plane wave basis set, in eV)
# PREC   =  Accurate   (Precision level: Normal or Accurate, set Accurate when perform structure lattice relaxation calculation)
LWAVE  = .TRUE.        (Write WAVECAR or not)
LCHARG = .TRUE.        (Write CHGCAR or not)
ADDGRID= .TRUE.        (Increase grid, helps GGA convergence)
LASPH  = .TRUE.        (Give more accurate total energies and band structure calculations)
PREC   = Accurate      (Accurate strictly avoids any aliasing or wrap around errors)
# LVTOT  = .TRUE.      (Write total electrostatic potential into LOCPOT or not)
# LVHAR  = .TRUE.      (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# LPLANE = .TRUE.      (Real space distribution, supercells)
# NWRITE = 2           (Medium-level output)
# KPAR   = 2           (Divides k-grid into separate groups)
# NGXF    = 300        (FFT grid mesh density for nice charge/potential plots)
# NGYF    = 300        (FFT grid mesh density for nice charge/potential plots)
# NGZF    = 300        (FFT grid mesh density for nice charge/potential plots)
#LDIPOL =  .TRUE.       (Dipole correction)
#IDIPOL =  3           (The direction for dipole correction)
#DIPOL  = 0.5 0.5 0.5 (The position for dipole correction)

Electronic Relaxation
ISMEAR =  1            (Gaussian smearing, metals:1)
SIGMA  =  0.1         (Smearing value in eV, metals:0.2)
NELM   =  999           (Max electronic SCF steps)
NELMIN =  6            (Min electronic SCF steps)
EDIFF  =  1E-05        (SCF energy convergence, in eV)
# GGA  =  PS           (PBEsol exchange-correlation)

Ionic Relaxation
NSW    =  999          (Max ionic steps)
IBRION =  2            (Algorithm: 0-MD, 1-Quasi-New, 2-CG)
ISIF   =  2            (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)
EDIFFG = -2E-03        (Ionic convergence, eV/AA)
# ISYM =  2            (Symmetry: 0=none, 2=GGA, 3=hybrids)

3. KPOINTS
0
Monkhorst-Pack
2    2    1
0.0  0.0  0.0
现在不太了解KPOINTS，这是文献设置值
这是我目前思考的程度，如果有理解错的地方也恳请看到帖子的朋友批评指正。

弛豫计算顺利完成，由于活性氧和金属表面作用较强，Li原子位置发生了一些改变，弛豫时间还是比较长的，主要是为了达到能量收敛标准。现在开始做频率计算。频率计算的四个输入文件准备如下、
1. POSCAR
将弛豫计算的CONTCAR直接拷贝变成POSCAR，注意删掉最下面的速度信息。
2. INCAR
处理做频率的相关算法设置，其余要与弛豫计算一致
Global Parameters
ISTART =  1            (Read existing wavefunction, if there)
ISPIN  =  2            (Spin polarised DFT)
NELECT =   167
NUPDOWN = 1
MAGMOM =  40*0.01 1.0 -1.0 3*0.01 3*0.01 4*0.01
ICHARG =  0            (Non-self-consistent: GGA/LDA band structures)
LREAL  = Auto       (Projection operators: automatic)
ENCUT  =  520          (Cut-off energy for plane wave basis set, in eV)
# PREC   =  Accurate   (Precision level: Normal or Accurate, set Accurate when perform structure lattice relaxat
LWAVE  = .TRUE.        (Write WAVECAR or not)
LCHARG = .TRUE.        (Write CHGCAR or not)
ADDGRID= .TRUE.        (Increase grid, helps GGA convergence)
LASPH  = .TRUE.        (Give more accurate total energies and band structure calculations)
PREC   = Accurate      (Accurate strictly avoids any aliasing or wrap around errors)
# LVTOT  = .TRUE.      (Write total electrostatic potential into LOCPOT or not)
# LVHAR  = .TRUE.      (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# LPLANE = .TRUE.      (Real space distribution, supercells)
# NWRITE = 2           (Medium-level output)
# KPAR   = 2           (Divides k-grid into separate groups)
# NGXF    = 300        (FFT grid mesh density for nice charge/potential plots)
# NGYF    = 300        (FFT grid mesh density for nice charge/potential plots)
# NGZF    = 300        (FFT grid mesh density for nice charge/potential plots)
#LDIPOL =  .TRUE.       (Dipole correction)
#IDIPOL =  3            (The direction for dipole correction)
#DIPOL  = 0.5 0.5 0.5 (The position for dipole correction)

Electronic
ISMEAR =  1
SIGMA  =  0.1         (please check the width of the smearing)
ALGO = Normal
EDIFF = 1E-5

Frequence Calculations
NSW    =  1            (number of ionic steps. Make it odd.)
IBRION =  5            (frequence calculation)
POTIM  =  0.02         (displacement step)
NFREE  =  2            (displacement freedom)
3. KPOINTS
设置为Gamma only的KPOINTS文件即可