2024-09-23

算H₂结合能

H₂结合能 = -4.53627867 eV

1. H原子

INCAR

System = H atom

! init
ISTART = 0; ICHARG = 2 

! large cell
PREC=Accurate

! electronic self-consistency
ENCUT = 300
ISMEAR = 0; SIGMA = 0.1
EDIFF = 1E-05

! spin polarization
ISPIN = 2

不太明白ISMEAR是什么，看了一篇官网的说明，里面提到：

The only difference to the bulk calculation is that Gaussian smearing must be used now. You might set SIGMA to a very small value; this is necessary if atomic orbitals are almost degenerated.

所以跟着设置了ISMEAR = 0; SIGMA = 0.1；ENCUT参照了POTCAR里的ENMAX值；有极性所以设置ISPIN = 2。

POSCAR

H atom
10
1 0 0
0 1 0
0 0 1
H
1
cart
0 0 0

KPOINTS

H atom
0
G
1  1  1
0  0  0
 

POTCAR

PAW_PBE H 15Jun2001                    
   1.00000000000000     
…… …… …… ……

得到结果，其中

  …… …… …… ……

  atomic energy  EATOM  =        12.44577538
  Solvation  Ediel_sol  =         0.00000000
  ---------------------------------------------------
  free energy    TOTEN  =        -1.11159491 eV

  energy without entropy =       -1.11159491  energy(sigma->0) =       -1.11159491

  …… …… …… ……

2. H₂分子

INCAR

System = H2 molecule

! init
ISTART = 0; ICHARG = 2 

! large cell
PREC=Accurate

! electronic self-consistency
ENCUT = 550
ISMEAR = 0; SIGMA = 0.1
EDIFF = 1E-5

! ionic relaxation
IBRION = 2
NSW = 50
EDIFFG = -0.01
ISIF = 2

看了官网这页的表格:

ISIF	calculate		degrees-of-freedom
	forces	stress tensor	positions	cell shape	cell volume
0	yes	no	yes	no	no
1	yes	trace only	yes	no	no
2	yes	yes	yes	no	no
3	yes	yes	yes	yes	yes
4	yes	yes	yes	yes	no
5	yes	yes	no	yes	no
6	yes	yes	no	yes	yes
7	yes	yes	no	no	yes
8	yes	yes	yes	no	yes

设置了ISIF = 2。

POSCAR

H2 molecule in a box
16
1 0 0
0 1 0
0 0 1
H
2
cart
0 0 0
0.1 0 0

KPOINTS

POTCAR

PAW_PBE H 15Jun2001                    
   1.00000000000000    

…… …… …… ……

得到的结果是

 FREE ENERGIE OF THE ION-ELECTRON SYSTEM (eV)
  ---------------------------------------------------
  free  energy   TOTEN  =        -6.77191945 eV

  energy  without entropy=       -6.77191945  energy(sigma->0) =       -6.77191945

CONTCAR

H2 molecule in a box                    
   16.0000000000000     
     1.0000000000000000    0.0000000000000000    0.0000000000000000
     0.0000000000000000    1.0000000000000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000    1.0000000000000000
   H 
     2
Direct
  0.0265643014703201 -0.0000000000000000  0.0000000000000000
  0.0734356985296813  0.0000000000000000 -0.0000000000000000

键长为 16 * (0.0734356985296813 - 0.0265643014703201) = 0.74994235295 ≈ 0.75

H₂结合能为 -6.77191945 - 2*（-1.11159491）= -4.54872963 eV

有关文献¹给出的参考值是-4.58eV

For all calculations, we use the generalized gradient corrections (GGA) of Perdew and Wang^31,32 commonly refered to as PW91. …… …… For the H2 dimer, we find a bond length of a0=0.750 A˚ , a binding energy of E=4.58 eV and a vibrational frequency of w=4300 cm−1.

3. 测试

测试了H原子的ENCUT和POSCAR里的比例系数；对H₂分子只测试了ENCUT。尝试写了脚本：

测试H原子ENCUT

#!/bin/bash
n=1
for i in $(seq 250 25 600) 650 700 750;
do
    cp -r H_atom/ H_1_$n/
    sed -i "/ENCUT/c\ENCUT = $i" H_1_$n/INCAR
    cd H_1_$n/
    qsub wz.pbs
    cd ..
    n=$((n+1))
done

测试H原子POSCAR

n=1
for i in $(seq 8 1 19);
do
cp -r H_atom/ H_2_$n/
cd H_2_$n/
cat > POSCAR <<HERE
H atom
$i
1 0 0
0 1 0
0 0 1
H
1
cart
0 0 0

HERE
qsub wz.pbs
cd ..
n=$((n+1))
done

测试H₂分子ENCUT

#!/bin/bash
n=1
for i in $(seq 250 25 600) 650 700 750;
do
    cp -r H2_molecule/ H2_molecule_$n/
    sed -i "/ENCUT/c\ENCUT = $i" H2_molecule_$n/INCAR
    cd H2_molecule_$n/
    qsub wz.pbs
    cd ..
    n=$((n+1))
done

测试结果在这里。取原子ENCUT = 475，盒子大小12，分子ENCUT = 500，得到结合能为 -6.77014457 -2 * (-1.11693295) = -4.53627867 eV

4. 问题

看PPT里算的是-4.497922eV，用里面的参数算了一遍，H原子得到了相同的结果，H₂分子走了一步就结束了，CONTCAR和POSCAR一模一样，得到的energy = -0.11878548(而PPT里却写算出了-6.731532eV，不太理解为什么会出现这种情况），如果像上面自己算的时候一样，额外设置EDIFF和EFDIFFG比默认值小一级，会得到-6.77041931（和之前算的结果差不多）
本来理解的是，direct要乘第二行的scailing factor，cartisian填写实际坐标，但看了官网里关于POSCAR的cartisian和direct的解释，好像是理解反了？

“Direct” means the positions are provided in direct (fractional) coordinates. “Cartesian” specifies that positions are provided in a Cartesian coordinate system. However, the actual ion positions are also multiplied with the universal scaling factor, i.e.

5. Teacher’s review

ISMEAR是对部分占据轨道的不同处理算法，因此区别主要体现在材料的输运性上，例如轨道部分占据常出现在金属材料中，因此在计算金属和非金属体系时需要设置不同的ISMEAR值，SIGMA可以看做是积分的能量展宽，在ismear取-5时，sigma为零。在对计算体系进行优化的过程中（不知道体系到底是金属还是半导体、绝缘体），一般多采用ISMEAR=0(Gaussian smearing)，配合设置一个较小的sigma(0.05)，即可。优化收敛后，需要计算对比total energy的时候，推荐采用ismear=-5, sigma=0

EDIFF和EDIFFG这两个参数一个作为体系优化收敛的能量判据，一个是收敛力判据(EDIFFG<0时)，两者要配合使用，不合理搭配会出现优化不易收敛的问题。提高EDIFF和EDIFFG的收敛精度，确实会导致最终收敛能量略低于之前的结果，但是会消耗时间降低效率，这个在实际操作中看情况选择即可。如果需要用到高精度收敛，建议分步优化，先相对低收敛精度优化，再提高优化收敛精度(同时调整优化IBRION–优化方式和POTIM–优化步幅/弛豫时长)

POSCAR的结构是：
第一行类似注释，比如这个结构是什么材料、组成、甚至比例，完全看你个人
第二行：scaling factor，比例系数，一般我们默认为1，对于需要整体膨胀或压缩的结构，在这里直接修改比例系数会比较方便；
第三行~五行：晶体笛卡尔坐标系；
第六行：POSCAR中元素
第七行：与元素对应的原子个数
第八行：Direct、Cartesian绝对
第九行以下的具体原子坐标采用的是分数坐标(Direct)还是笛卡尔坐标也就是绝对坐标(Cartesian)。
这次你测试H2用的立方晶格比较特殊，可以直接通过scaling factor来同时修改晶格常数，后续实际材料晶格常数未必相同，建议后续scaling factor默认取1即可。

*例如，这两个等效：（1）第二行1，第三行800，第八行Direct，第九行的x是0.25；（2）第二行1，第三行800，第八行Cartesian，第九行的x是2；

G. Kresse and J. Hafner, Surface Science 459, 287 (2000). ↩

算H2结合能

ToC