| Martin's profileMartin's FactoryPhotosBlogLists | Help |
|
September 26 赝势方法;模型势方法;pseudo-potential method 分子式: CAS号: 性质:又称模型势方法。即价电子从头计算法。其基本思想是,所研究体系的哈密顿算符仅显含价电子部分,而将原子内层的全部电子连同原子核构成的核实对外层价电子的作用用适当的模型势函数(即赝势pseudo potential)表示,同时引入表示投影算符的势以便将价电子波函数与内层电子波函数分离开来,然后对价电子进行变分并用自洽迭代处理,计算出价电子轨道波函数和能级值等。所引入的两个势均通过全电子的原子从头计算确定,随后用于分子计算。常见的赝势有两类,一类是基于固体物理理论的菲利普斯-克兰曼(Phillips-Kleinman)赝势原理发展起来的,一类是从早期的价电子概念出发,由赫律纳加(Huzinaga)等提出的模型。赝势方法主要用于含有重元素的化合物或原子簇的研究。 VASP赝势文件POTCAR详解可以理解为分子力学模拟中的力场文件 但包括的信息更多 VASP4.6将各元素优化的INCAR里的参数也包括在这里了,作为支持PREC的缺省选择 通常各元素的POTCAR已经包括在软件包里了 我们只需要按照POSCAR里的顺序,将各元素的POTCAR按顺序连接起来就可以了 如以下命令: cat file1 file2 file3 > POTCAR 软件包自带的绝大多数赝势是超软赝势(US-PP)了,但不少元素有两个版本,如何 选取呢? 一个简单的办法是看后缀 标准的没有后缀 _h 硬一点 _s 软一点 _pv,_sv,_d 就是说semi-core的p,s或者d也当做价态处理了 如果是数字的话,表示的可能是不同的半径截距 也可以参考各版本同目录下的V_RHFIN file ,PSCTR file 这两个文件告知该版本的赝势是如何生成的。比如: V_RHFIN file Sc: 6p d2 s1 8 21. .002000 44.95590 125. .25E-05 .300 200FCA 12.00000 .7 1.0 0 1.0 .0 .5 -320.8847 2.0000 2.0 .0 .5 -34.4217 2.0000 2.0 1.0 1.5 -28.2366 6.0000 3.0 .0 .5 -3.7944 2.0000 3.0 1.0 1.5 -2.2591 6.0000 3.0 2.0 2.5 -.1113 2.0000 4.0 .0 .5 -.2699 1.0000 4.0 3.0 2.5 -.1000 .0000 第一行是注释行 给出基本的信息 第二行是最重要的控制行 8 21. .00 2000 44.95590 125. .25E-05 .300 200 F CA 12.00000 J Z XION N AM H DELRVR PHI NC1 | CH QCOR | GREEN J - 轨道数 Z - 原子序数 XION - 离子化程度 一般设为0 N - 格点数 AM - 原子质量 H - 决定格点间距 DELRVR - 自洽收敛标准 PHI - 线性拟合参数 NC1 - 最大自洽循环次数 GREEN - 是否存在初始的势 CH - 交换相关能(XC)类型 Slater-XC HL Hedin Lundquist (1971) CA Ceperly and Alder parameterized by J.Perdew and Zunger WI Wigner interpolation PB Perdew -Becke PW Perdew -Wang 86 LM Langreth-Mehl-Hu 91 Perdew -Wang 91 QCOR - 非价键电子数(core electrons) 第三行开始是每个轨道的具体参数,依次为 n l j(=l±1/2) 原子轨道能 占有率 PSCTR file of LDA/H1.25 TITEL = US H LULTRA = T use ultrasoft PP ? RWIGS = 0.57 nn distance ! Wigner-Seitz radius RCLOC = .65 NE = 100 LCOR = .TRUE. QCUT = -1 RMAX = 3.0 ! core radius for proj-oper Description l E TYP RCUT TYP RCUT(cutoff radius) 0 0 15 0.80 23 1.25 0 0.5 15 0.80 23 1.25 1 -0.2 15 0.80 23 1.25 最重要的地方上面已经用颜色标出来啦:) 说明一下,TYP是指赝势的类型,RCUT是半径截距,TYP可取的值如下: 正则 1 BHS 2 TM 3 VAN 6 XNC 7 RRKJ wave function possibly with node 15 RRKJ wave function strictly no node 非正则 +8 最后一个问题是LDA or GGA。貌似没有定论目前。 这个最好是两个一起做做看啦。或者看文献别人验证过哪个数据好。 其实据说目前最好的是PAW(P.E.Blochl,Phys.Rev.B 50,17953(1994).,Phys.Rev.B 59,1758(1999).)。 VASP新手入门十个简单问题1.TEBEG,TEEND的单位是K吗? 答:是K为单位的。 2.在运算中,直接COPY POSCAR文件到工作目录就可以吗,用不用再POSCAR或INCAR文件中设定参数来控制POSCAR的赝势的选取方法? 答:不大明白你的问题。在一个工作目录准备好INCAR,POSCAR,POTCAR,KPOINTS这四个文件。先看POTCAR中LEXCH后面的值是什么,如果是CA,那么INCAR中就不用设置GGA了;如果是91,那么就在 INCAR 中输入GGA=91;如果是PE,那么就设置GGA=PE。 3.模拟表面,一般SLAB选多厚? 答:自己多测试不同的厚度后,来确定slab 的厚度 4.在KPOINT文件中,如果选取MONKHORST-pack的K点产生办法,K点相对自动网格移动时怎么回事? 答:你可以看看monkhorst,pack1976年在phys. rev. b上的文献,他们在文献中给出了对于sc,bcc和fcc三种格子的一种效率很高的k点选取标准。但是这种k-mesh的选取方法并不是将第一个点取在倒空间的原点处,而是有一个shift,比如对于fcc而言,这个shift时(1/8,1/8, 1/8)。但是这种mp方法对于正交格子适用,而对于三方或者六方的格子则并不好。因此这时建议采用\gamma点为中心的k-mesh(vasp用户手册) 5.想做表面重构,VASP 能做吗,能用分子动力学模拟吗,怎么设置参数? 答:做重构的计算,一般用分子动力学的方式来模拟。你可以参考文献上模拟表面体系重构的作法。大多是采用比较几个表面模式的能量以及相关的热力学量。 6.不指定电子的总SPIN,会有影响吗,一般 什么体系才考虑自旋? 答:你做的体系是否包含3d-过渡金属离子(V至Ni)或其他含f电子的原子,如果没有的话,可以不考虑自旋极化的计算。 7.表面的原子或分子的运动,扩散还有重构,IBRION怎么设定? 答:是扩散的化,一般用NEB的方法计算其扩散路径或扩散势垒。 8.如果知道一个化合物的组成,想计算它的总能一般应怎么办?同时如果一直组成如何在常温下计算它的稳定结构,比如FeMgCl2? 答:计算总能,自己准备好vasp的输入文件,运算后就直接得到。可能考虑它的结合能比较好。 9.在计算一个系统的稳定结构时,SMASS参数如何设置? 答:看你是计算这个稳定结构的什么性质:计算总能和DOS,用ISMEAR=-5,不用管SMASS。 10.具体给介绍一下STM图像的模拟办法,如何计算特定的能带,还有计算完成后一般去那个文件查找对自己有用的数据? 答:一般来说用parchg做个图就可以了,当然以上方法仅限于T-H近似。 如何用VASP计算晶格常数我们用Pd金属作为例子。 如何用VASP计算单个原子的能量和能级 氢原子的能量为-13.6eV在这一节中,我们用VASP计算H原子的能量。 对于原子计算,我们可以采用如下的INCAR文件 PREC=ACCURATE NELMDL = 5 make five delays till charge mixing ISMEAR = 0; SIGMA=0.05 use smearing method 采用如下的KPOINTS文件。由于增加K点的数目只能改进描述原子间的相互作用,而在单原子计算中并不需要。所以我们只需要一个K点。 Monkhorst Pack 0 Monkhorst Pack 1 1 1 0 0 0 采用如下的POSCAR文件 atom 1 15.00000 .00000 .00000 .00000 15.00000 .00000 .00000 .00000 15.00000 1 cart 0 0 0 采用标准的H的POTCAR 得到结果如下: k-point 1 : 0.0000 0.0000 0.0000 band No. band energies occupation 1 -6.3145 1.00000 2 -0.0527 0.00000 3 0.4829 0.00000 4 0.4829 0.00000 我们可以看到,电子的能级不为-13.6eV。 Free energy of the ion-electron system (eV) --------------------------------------------------- alpha Z PSCENC = 0.00060791 Ewald energy TEWEN = -1.36188267 -1/2 Hartree DENC = -6.27429270 -V(xc)+E(xc) XCENC = 1.90099128 PAW double counting = 0.00000000 0.00000000 entropy T*S EENTRO = -0.02820948 eigenvalues EBANDS = -6.31447362 atomic energy EATOM = 12.04670449 --------------------------------------------------- free energy TOTEN = -0.03055478 eV energy without entropy = -0.00234530 energy(sigma->0) = -0.01645004 我们可以看到TOTEN-EATOM也不等于-13.6eV。 在上面的计算中有个问题,就是H原子有spin,而在上面的计算中我们并没有考虑到spin。 所以如果我们改用LSDA近似,在INCAR中用ISPIN=2的tag,则得到如下结果: k-point 1 : 0.0000 0.0000 0.0000 band No. band energies occupation 1 -7.2736 1.00000 2 -0.1229 0.00000 3 0.4562 0.00000 4 0.4562 0.00000 5 0.4562 0.00000 spin component 2 k-point 1 : 0.0000 0.0000 0.0000 band No. band energies occupation 1 -2.4140 0.00000 2 -0.0701 0.00000 3 0.5179 0.00000 4 0.5179 0.00000 5 0.5179 0.00000 Free energy of the ion-electron system (eV) --------------------------------------------------- alpha Z PSCENC = 0.00060791 Ewald energy TEWEN = -1.36188267 -1/2 Hartree DENC = -6.68322940 -V(xc)+E(xc) XCENC = 2.38615430 PAW double counting = 0.00000000 0.00000000 entropy T*S EENTRO = 0.00000000 eigenvalues EBANDS = -7.27361676 atomic energy EATOM = 12.04670449 --------------------------------------------------- free energy TOTEN = -0.88526212 eV energy without entropy = -0.88526212 energy(sigma->0) = -0.88526212 氢原子的能量约等于-12.92eV。可以看到在LDA中如果限制自旋,使能级大概提高了0.88eV。 如果我们采用GGA的赝势,并且同样打开自旋限制(ISPIN=2),在此例子中,得到的结果将更加精确 Free energy of the ion-electron system (eV) --------------------------------------------------- alpha Z PSCENC = 0.00621465 Ewald energy TEWEN = -1.38027565 -1/2 Hartree DENC = -6.91107031 -V(xc)+E(xc) XCENC = 2.32601856 PAW double counting = 0.70286470 -0.71256934 entropy T*S EENTRO = 0.00434843 eigenvalues EBANDS = -7.71257941 atomic energy EATOM = 12.52153358 --------------------------------------------------- free energy TOTEN = -1.15551480 eV energy without entropy = -1.15986323 energy(sigma->0) = -1.15696428 用GGA算得的氢原子的能量约等于-13.67eV。 但是如何理解所得到的能级,由于用到了赝势,本人并不很清楚如何解释能级意义,欢迎大家指教。 如何用VASP计算单个原子的能量[权威版] 这是VASP作者给出的解答,所以算是权威了.原文如下: The energy zero .. in VASP all energies are referred to the the reference state for which the potential was generated! this is in most cases not the real groundstate of the atom .. to determined the energy of the grounstate of the atom place the atom in a larger non cubic box to break initial symmetry (i.e. 11Ax10Ax9A) use the G point only INCAR: ISPIN = 2 ! spin polarized ISMEAR = 0 ; SIGMA = 0.2 ! for small sigma conv. for TM is diff. MAGMOM = 2 ! initial magnetic moment one should use the energy value energy without entropy of the OUTCAR file since this coverges most rapidly to the correct energy for sigma 0 .. Ecoh=E_total-n*Eatom 大约翻译如下: 在VASP中,所有的能量是相对于产生赝势的组态来说的,所以并不能将计算得到的能量当成真正的原子的基态能量。 为了得到原子的基态能量,将原子放在一个非正方体的晶格内,比如说11Ax10Ax9A的晶格(这样做是为了消除简并)。 KPOINTS只用G撒点 大概的INCAR文件可以如下: ISPIN = 2 ! spin polarized ISMEAR = 0 ; SIGMA = 0.2 ! for small sigma conv. for TM is diff. MAGMOM = 2 ! initial magnetic moment 计算结束后,应该用OUTCAR中不含entropy的那一项能量,因为它通常更快的收敛到sigma=0的正确值。最后计算结合能的公式应该为 Ecoh=E_total-n*Eatom 用VASP计算氢分子的能量 氢气分子的解离能,也就是结合能,根据资料中给出的是约4.48eV。(G. Kresse & J. Hafner, Surface Sci. 459 (2000) 287) 为此,首先要计算一个氢气分子的孤立能量,再减去两个孤立氢原子能量,将得到氢气分子的结合能。 先计算单个原子能量,选取PAW_PBE文件夹下的H下面的POTCAR 用到的其他输入文件如下: INCAR: SYSTEM = H atom in a box ISMEAR = 0 ! Gaussian smearing SIGMA = 0.01 ENCUT = 350.0 KPOINTS: Automatic mesh 0 Monkhorst Pack 1 1 1 0. 0. 0. POSCAR: H atom in a box 1.0 ! universal scaling parameters 7.0 0.0 0.0 ! lattice vector a(1) 0.0 8.0 0.0 ! lattice vector a(2) 0.0 0.0 9.0 ! lattice vector a(3) 1 ! number of atoms cart ! positions in cartesian coordinates 0 0 0 用上述文件计算得到 TOTEN(H) = +0.000854eV (H atom),用PAW_GGA得到类似结果。 然后计算氢分子能量,用类似的输入文件: INCAR: SYSTEM = H2 dimer in a box ISMEAR = 0 ! Gaussian smearing NSW = 5 ! 5 ionic steps IBRION = 2 ! use the conjugate gradient algorithm ENCUT = 350.0 POTIM = 0.1 KPOINTS: Automatic mesh 0 Monkhorst Pack 1 1 1 0. 0. 0. POSCAR: H2 molecule in a box 1.0 ! universal scaling parameters 8.0 0.0 0.0 ! lattice vector a(1) 0.0 8.0 0.0 ! lattice vector a(2) 0.0 0.0 8.0 ! lattice vector a(3) 2 ! number of atoms cart ! positions in cartesian coordinates 0 0 0 ! first atom 0 0 0.5 ! second atom 根据定义,E=-(TOTEN(H2)-2*TOTEN(H)) = 6.68eV,键长 = 0.75102 A。比文献中的结合能大了不少。 出现上述结果的原因是在计算单个H原子能量的时候没有指定基态为spin polarized state。为了得到正确的解离能,计算单个H原子能量的时候需要指定ISPIN=2。加上以上的tag后,TOTEN(H atom [spin-polarized])=-1.10351 eV,用公式重新计算解离能: E_binding=-(TOTEN(H2)-2*TOTEN(H [spin-polarized ])) = 4.48eV 和文献吻合。另外,文献[G.Kresse, PRB 62, 8295 (2000)]中详细讨论了氢分子的解离能的计算方法和结果。 用VASP进行Partial Charge分析实例VASP Version : 4.6 在这篇文章中,我将首先介绍Partial Charge的概念,以及如何用VASP具体的计算Partial Charge。首先,所谓的Partial Charge是针对与Total Charge来说的,指的是某个能量范围、某个K点或者某个特定的态所对应的电荷密度。在文献中最常见的是价带顶部,导带底部,表面态或者局域态所对应的Partial Charge。通过分析这些态所对应的Partial Charge,可以得到体系的一些性质,比如局域态具体的是局域在哪个原子上等。我将通过具体的例子说明如何用VASP进行Partial Charge Analysis。 进行Partial Charge Analysis的第一步是进行自洽的计算,得到体系的电子结构。这一步的计算采用通常的INCAR和KPOINTS文件。在自洽计算结束后,我们需要保存WAVECAR文件。(通过在INCAR文件中设置LWAVE=TRUE实现)在这个例子中,假设我们需要计算一个硅纳米线的导带和价带的Partial Charge。硅纳米线的结构如下:
http://www.quantumchemistry.net/Experience/UploadFiles_6872/200603/20060331155154521.jpg
http://www.quantumchemistry.net/Experience/UploadFiles_6872/200603/20060331155540648.jpg 第一种Partial Charge分析的INCAR LPARD=.TRUE. 这样的INCAR给出的是指定能带,指定K点所对应的Partial Charge。分析导带、价带等的Partial Charge特性,通常采用的都是这种模式。 第二种Partial Charge分析的INCAR LPARD=.TRUE. 这样的INCAR给出的是在[-10.3 -5.1]能量之间的Partial Charge。这种模式适合于分析某个能量区间内的波函数的性质。 第三种Partial Charge分析的INCAR LPARD=.TRUE. 这样的INCAR给出的是从[Ef-1.0 Ef]能量之间的Partial Charge。这种模式最利于分析费米面附近的波函数的性质。 用第一种方法,我们可以得到硅纳米线价带顶部和导带底部的Partial Charge如下:
http://www.quantumchemistry.net/Experience/UploadFiles_6872/200603/20060331155753136.jpg
|
|
|