第一性原理简介.doc

简介：ABINIT的主程序使用赝势和平面波，用密度泛函理论计算总能量，电荷密度，分子和周期性固体的电子结构，进行几何优化和分子动力学模拟，用TDDFT（对分子）或GW近似（多体微扰理论）计算激发态。此外还提供了大量的工具程序。程序的基组库包括了元素周期表1-109号所有元素。ABINIT适于固体物理，材料科学，化学和材料工程的研究，包括固体，分子，材料的表面，以及界面，如导体、半导体、绝缘体和金属。

功能：
可以计算很多物理属性：
A. 计算倒格子中核与电子的总能量。
A.1. 计算使用平面波和赝势。
A.2. 总能量的计算使用密度泛函理论（DFT）。可以使用大多数重要的局域密度近似 (LDA)，包括Perdew-Zunger近似。可以使用两种不同的局域自旋密度(LSD)，包括Perdew Wang 92和M. Teter的LSD。还可以使用Perdew-Burke-Ernzerhof，revPBE，RPBE和HCTH等GGA (自旋极化和非极化)。
A.3. 自恰场计算生成DFT基态，以及相关的能量和密度。此后的非自恰计算可以对能带结构的大量k-点产生本征能量。态密度的计算即可以用四面体方法，也可以用模糊技术。
A.4. 程序可以使用多种不同的赝势。对整个周期表适用的有两种：Troullier-Martins型和Goedecker型（这种类型包括自旋-轨道耦合）如果需要的话，有四个代码可以产生新的赝势。
A.5. 程序本身可以处理金属和绝缘体系。
A.6. 晶胞可以是正交或者非正交。计算可以输入任何对称性及相应的k-点集。
A.7. 电子体系可以用自旋极化和自旋非极化计算。一个特殊的选项可以有效地处理反铁磁性。可以对总能量计算非共线的磁性（不能用于力，张量，相应函数...）。可以禁止晶胞的总磁矩。
A.8. 总能量，力，张量和电子结构的计算可以考虑自旋-轨道耦合。
A.9. 能量可分解为不同的成分（局域势，XC，Hartree...）。
A.10. 计算内部电子本征值。
A.11. 230个空间群和1191个Shubnikov磁群的对称性分析。
B. 计算总能量和本征能量
B.1. 用解析公式计算Hellman-Feynman力。
B.2. 计算应力。
B.3. 极化的计算。
B.4. 响应的计算。
B.5. 计算近似的和准确的磁化系数矩阵和介电矩阵。
B.6. 解析计算电子本征能量的导数。
B.7. 计算光学传导性。
B.8. Born有效电荷的能带分解，以及局域化张量的计算。
C. 激发态
C.1. 用GW近似计算电离能和亲和能。
C.2. 用TDDFT计算原子和分子的（单重、三重）激发态和振荡强度。
D. 移动原子，改变晶胞参数
D.1. 用不同的方法寻找平衡构型。可以同时优化晶胞参数。优化过程中如果需要的话，可以固定指定的晶胞参数，角度，或原子位置。
D.2. 有两种算法进行分子动力学计算。
D.3. 自动分析键长键角。原子坐标的格式支持用可视化软件XMOL显示。
E. 分析和图形工具
E.1. 后期处理程序cut3d用于分析密度和势文件。它还可以改变文件格式，提取2D明面或者1D线。此外还可以分析波函文件。
E.2. 另一个后期处理程序aim，用于进行Bader的“原子中的分子”（AIM）密度分析。
E.3. 对可视化程序产生格式化数据：键结构（用XMGR显示），不同参数的总能量（用XMGR显示），电荷密度（3D轮廓线，先用cut3d，再用商业程序matlab；cut3d也可以产生2D密度图）。
E.4. 后期处理程序band2eps自动画出eps格式的声子散射曲线。

平台：各类Unix/Linux，Win98/Dos

相关软件：

1. FHI98PP

2. JMol

3. Xtal

4. XtalEdit

5. MAPS

6. http://theory.lsi.polytechnique.fr/codes/codes.html
一些GW程序，有的包含在ABINIT中，有的具有到ABINIT和PWSCF的接口。

ABINIT教程，第五课：AlAs的动力学和介电性质,  在这一课中，将介绍如何计算得到绝缘体材料的以下物理性质：
  在Γ点的声子频率和本征矢
  介电常数
  Born有效电荷
  LO-TO分裂
  Brillouin区非Γ点的声子频率和本征矢
  原子间力常数（还未完成，有待补充）
  从原子间力常数计算声子色散曲线（还未完成，有待补充）
  有关的热力学量（还未完成，有待补充）

  我们将学习使用ABINIT能计算响应函数的特点。在将来的版本中，我们将继续学习相关的代码Mrgddb和Anaddb部分。

  这一课需要花几个小时才能完成。

ABINIT计算实例分析2: 计算电子态密度

电子态密度(Density of States, DOS)定义为单位能量范围内所允许的能级数。在电子结构计算程序中进行计算时，一般是按如下的步骤来进行的：先采用小的k点网格进行自洽计算，然后采用较密的k点网格进行非自洽计算。

下面介绍采用ABINIT程序电子态密度时的步骤，解释需要输入的关键词，以及计算结果的处理说明。这里以采用ABINIT5.2.4版本计算MgO电子态密度（包括总态密度和spd分波态密度等)为例子。

在例子中MgO的晶格常数(它是面心立方晶体，NaCl结构，a=7.921 a.u.)[1]取为实验值，平面波切断动能和k点网格大小分别取为24 Ha，10x10x10（自洽计算中）和16x16x16（态密度计算中）。
[1] O.L. Anderson and P. Andreatch, J. Am. Ceram. Soc. 49, 404 (1966).

一、计算总态密度
通过abinit中的ndtset设置两组数据进行计算，先是自洽计算，后是态密度计算的设置。在计算态密度时，需将prtdos设置为1、2或3。当prtdos设置1时，采用离散(smearing)的方法来确定态密度，即由本征值、k点权重(参数wtk的值)、每个本征值的占有数（由occopt设置的参数所选择的方法来确定占有数)以及离散展宽参数(tsmear参数设置的值)来计算，此时需要注意或设置wtk、occopt和tsmear参数。一般对一维体系或者在k点较少的情况采用prtdos=1的方法来计算态密度，以得到较为光滑的DOS曲线。当设置prtdos=2时，表示采用四面体方法来计算总态密度，此时不需要设置occopt和tsmear了，所采用的k点需要由参数ngkpt或kptrlatt参数来确定。当prtdos=3时，表示采用四面体方法计算分波态密度(angular-momentum projected)，同时需设置原子球半径。

以计算面心立方MgO的总态密度为例子。
i). in.files输入文件的内容:
tdos.in
tdos.out
mgoi
mgoo
mgo
12mg.pspnc
8o.pspnc
ii)主要输入文件tdos.in的内容：
ndtset 2
kptopt 1
nshiftk1 4
shiftk1
0.5 0.5 0.5
0.5 0.0 0.0
0.0 0.5 0.0
0.0 0.0 0.5
ngkpt1 8 8 8

prtden1 1
toldfe1 1.0d-6
iscf1 5

#Dataset 2 : DOS calculation
iscf2 -3
getden2 1
getwfk2 1
prtdos2 2
dosdeltae 000005
ngkpt2 16 16 16
tolwfr2 1.0d-16
shiftk2
0.0 0.0 0.0
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0

#Definition of the unit cell
acell 3*7.921
rprim
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0

ntypat 2
znucl 12 8
natom 2
typat 1 2
xred
0.0 0.0 0.0
0.5 0.5 0.5

nband 12

ixc 1

ecut 24.0

nstep 45
diemac 9.0

这里利用ndtset设置两组数据，第一组是进行自洽计算以得到电荷密度，第二组是非自洽计算以得到总态密度。在第一组数据中:
kptopt 1 #设置通过ngkpt由程序自动生成k点
nshiftk1 4 #由于是fcc结构，为了使得的产生的k点对计算更有效，这里将所产生的k点进行移动
shiftk1 #设置了移动的量
0.5 0.5 0.5
0.5 0.0 0.0
0.0 0.5 0.0
0.0 0.0 0.5
ngkpt1 8 8 8 #设置k点网格的大小。由于是在自洽计算中，此时设置的数较小。

prtden1 1 #表示要将电荷密度输出到文件中，以供第二组中的计算用到。
toldfe 1.0d-6 #自洽迭代计算中总能变化收敛的标准
iscf1 5 #表示采用CG方法来优化总能以得到基态

在第二组数据中：
iscf2 -3 #表示进行非自洽计算，一般用在态密度和STM计算以及线性响应计算中对ddk微扰的计算。
getden2 1 #表示从上一步中得到的电荷密度作为初始的电荷密度读入
getwfk2 1 #表示读入上一步的波函数
prtdos2 2 #表示要输出总态密度（或进行总态密度计算)
dosdeltae 0.00005 #态密度计算时能量轴上的能量刻度
ngkpt2 16 16 16 #在态密度计算时，k点网格应该设置的密一些。
tolwfr2 1.0d-16
shiftk2
0.0 0.0 0.0
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0

其他的参数是用来设置晶格结构以及原子的种类和坐标，平面波切断动能，自洽收敛参数。

在计算出来的_DOS文件中给出了总态密度的值，其格式如下：
ABINIT package : DOS file

nsppol = 1, nkpt = 145, nband(1)= 12
Tetrahedron method
For identification : eigen(1:3)= -0.429 0.203 0.203

Fermi energy : 0.20340927

The DOS (in electrons/Hartree/cell) and integrated DOS (in electrons/cell) are computed,
at 38001 energies (in Hartree) covering the interval
between -0.5000 and 1.4000 Hartree by steps of 0.00005 Hartree.

index energy(Ha) DOS integrated DOS
0 -0.50000 0.0000 0.0000
1 -0.49995 0.0000 0.0000
2 -0.49990 0.0000 0.0000
..........
index给出的是能量点数的序号，第二列给出的是能量的数值（单位是Hatree），在画图是一般以eV，并减去前面给出的Fermi energy。第三列是总态密度值，第四列是电子数。

二、计算各个原子的分波态密度
ndtset 2
kptopt 1
nshiftk1 4
shiftk1
0.5 0.5 0.5
0.5 0.0 0.0
0.0 0.5 0.0
0.0 0.0 0.5
ngkpt1 8 8 8

prtden1 1
toldfe1 1.0d-6
iscf1 5

#Dataset 2 : DOS calculation
iscf2 -3
getden2 1
getwfk2 1
prtdos2 3
dosdeltae 0.00005
natsph2 2
ratsph2 2.880 2.880
ngkpt2 16 16 16
tolwfr2 1.0d-16
shiftk2
0.0 0.0 0.0
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0

#Definition of the unit cell
acell 3*7.921
rprim
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0

ntypat 2
znucl 12 8
natom 2
typat 1 2
xred
0.0 0.0 0.0
0.5 0.5 0.5

nband 12

ixc 1

ecut 24.0

nstep 45
diemac 9.0

在进行分波态密度计算时，与总态密度计算不同的是：prtdos须设置为3，另外需指定下面两个参数：
natsph 表示要计算多少个原子的分波态密度。最好设置为原子的总个数，实际上它的默认值是natom（即原子的总个数）。
ratsph 设置每个原子的态密度积分时的半径。由于现在只能对所有原子的ratsph都指定为一样的数，那么要计算不同原子的分波态密度，要设置为不同的ratsph值再计算一次。

  第5课的内容

  5.1 计算AlAs的基态几何结构和性质
  5.2 采用冷冻声子的方法计算总能的二阶导数
  5.3 采用响应函数的方法计算总能的二阶导数
  5.4 采用响应函数的方法计算在Γ点的动力学矩阵
  5.5 采用响应函数的方法计算体系对均匀电场的响应
  5.6 采用响应函数的方法计算非零q点的声子频率
  5.7 计算整个声子色散曲线以及相关的热力学量

  5.1计算AlAs的基态几何结构和性质

  先创建一个新的工作目录，名为work5，在该目录下准备输入文件：t5.files，t51.in以及Al和As的赝势文件。其中t5.files中列出了计算所用到的一些输入文件以及包含了计算结果的输出文件的文件名；t51.in文件中的内容是计算中最重要的一个输入文件，它控制了计算的流程和内容，定义了计算中各种输入参数的赋值。在这里列出了t5.files文件中的内容：

  ----t5.files的内容---- ----注释---

  t51.in #主要的输入文件是t51.in，定义了关键词的赋值

  t51.out #把主要的计算结果写在该文件中

  t51i #指明输入文件是以t51作为文件名开头的

  t51o #指明输出文件是以t51作为文件名开头的

  t51 #指明计算中有关的输入文件是t51作为开头的

  ../Psp/13al.981214.fhi #指明Al的赝势文件所在的位置（相对路径）

  ../Psp/33as.pspnc #指明As的赝势文件所在的位置

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

  同样也列出了t51.in文件中的内容：

  ----t51.in的内容-------

  # 计算晶体AlAs的总能

  #

  #指定进行基态计算

   kptopt 1 # 程序自动产生k点，并考虑对称性

   tolvrs 1.0d-18 # 自洽循环计算停止的标准

   iscf 5 # 表示采用算法5进行自洽循环计算

  #######################################################################

  #普通的输入参数

  #定义原胞

  acell 3*10.61 # 这等价于定义10.61 10.61 10.61

  rprim 0.0 0.5 0.5 # 这些同acell一起定义了基矢 ，无默认值

   0.5 0.0 0.5

   0.5 0.5 0.0

  #定义原胞中原子的类型

  ntypat 2 # 表示原胞中有2类原子

  znucl 13 33 # 关键词znucl用来指明每类原子的原子序数，其顺序必须与在上

   面的’files’文件中列出的赝势文件相一致，这里表明第一类原子

   是Al，第二类是As

  #定义原胞中的原子

  natom 2 # 表明原胞中有2个原子

  typat 1 2 # 第一个原子是属于第一类的(Al)，第二个原子是输入第二类的(As)

  xred # 这个关键词按下面的格式用来指明原子的位置，每个原子有三个数

   0.0 0.0 0.0 # 这三个数给出了第1个原子的分数坐标

   0.25 0.25 0.25 #这三个数给出了第2个原子的分数坐标

  #明确给出所要计算的能带数（最好不采用默认值）

  nband 4 # 对绝缘体（且要是DFT能正确计算出真正是绝缘体），则在响应

   函数的计算中没有必要考虑计算导带

  #定义交换关联函数

  ixc 1 # 采用Teter Pade参数化的LDA

  #定义平面波基矢

  ecut 3.0 # 平面波切断动能为3.0 Hartree

  #定义k点网格的设置

  kptrlatt -4 4 4 # 按Cartesian坐标来设置k点网格，这个网格是立方的而且它对

   4 -4 4 # 应的是所谓的8 x 8 x 8 Monkhorst-Pack 网格。当然也可以通过采

   4 4 -4 # 用ngkpt, nshiftk和shiftk来设置

  #定义SCF计算中的参数

  nstep 15 # SCF循环的最大步数

  diemac 9.0 # 尽管这个不是必须的，对预处理SCF循环是很值得的。模型介点

   函数被用标准的预处理，见dielng输入参数的描述。AlAs的介点

   常数要比Si（＝12）的小。