xiaoqiugood的个人博客分享 http://blog.sciencenet.cn/u/xiaoqiugood

博文

VASP系列:vdW、HSE、magenetic、Raman谱之vdW计算

已有 34200 次阅读 2013-10-11 02:35 |个人分类:关注的问题|系统分类:科研笔记

关注问题:

1) 什么情形下,需考虑分子间相互作用vdW?  vdW的计算原理及过程;

2) 能带结构精确计算的HSE方法原理及参数设置

3) 什么情形下需考虑magenetic?  magenetic计算的原理及参数设置;

4)  什么是自旋轨道耦合,什么情形下需考虑该作用;参数如何设置;

5)  Ramman谱的计算原理及振动分析方法


题记:Thanks must be given to Andreas

   做学问与做人: 我们不仅仅要学习导师的理论知识,还应该学习怎样做人;学习怎样雪中送炭,怎样诲人不倦。

   其实,有时懂与不懂、会与不会之间隔着的就是一层纸,如果有人点破则可以少走很多弯路,省很多时间。 可是,不是人人都会诲人不倦的,不是人人都具有送人玫瑰、手留余香的精神。

  有的人不仅不给你点破,而且还会故意在你和真理之间设置层层面纱,这或许就是人品的问题。中国有很多Nice的课题组,这些课题组都有很Nice的自己开发的软件,但有几个可以做到像calypso课题组一样?!

  Mulitifuntion的开发者sobereva则是另一个楷模。


1.来自Andreas的信件


You attached theScheffler vdW-paper, and it is quite right: the transition pressures betweenthe low-pressure molecular ice phases are much better described if vdW-effects are included in thecalculation. I found that of the VASP-supported correction schemes, the vdw-DF approach(implemented by Klimes) is very good to reproduce experimental ground statestability ranges of ice phases.

I have not tested this forH2O2. For the technical side, I have attached an INCAR file to use this schems(needs VASP 5.1+), INCAR.relax_vdw.

I have also attached two INCARfiles for a HSE hybrid-DFT calculation, INCAR.energy and INCAR.relax_hse.

Use INCAR.energy first to create DFT wave functions, then use themas starting guess for the HSE part. The VASP manual onlineexplains this procedure as well.


(1) INCAR.relax_vdw



# output options

LWAVE  = .FALSE. # write or don't write WAVECAR

LCHARG = .FALSE. # write or don't write CHG and CHGCAR

LELF   = .FALSE. # write ELF


# ionic relaxation

NSW = 100        # number of ionic steps

IBRION = 1       # 2=conjucate gradient, 1=Newton like

ISIF = 3         # 3=relax everything, 2=relax ions only, 4=keep volume fixed


# precision parameters

EDIFF = 1E-7     # 1E-3 very low precision for pre-relaxation, use 1E-5 next

EDIFFG = -1E-3    # usually: 10 * EDIFF

PREC = high      # precision low, med, high, accurate


# electronic relaxation

ISMEAR = 0       # -5 = tetraedon, 1..N = Methfessel

SIGMA = 0.1

ENCUT = 800      # cutoff energy

PSTRESS = 0


# Dion/Soler/Klimes vdW-DF

LUSE_VDW = .TRUE.

AGGAC = 0.0000

GGA = BO

PARAM1 = 0.1833333333

PARAM2 = 0.2200000000




2.vdW-DF functional of Langreth and Lundqvist et al.

http://cms.mpi.univie.ac.at/vasp/vasp/vdW_DF_functional_Langreth_Lundqvist_et_al.html


The vdW-DF proposed by Dion et al. [122] is a non-local correlation functional that approximately accounts for dispersion interactions. In VASP the method is implemented using the algorithm of Roman-Perez and Soler [123] which transforms the double real space integral to reciprocal space and reduces the computational effort. Several propsed versions of the method can be used: the original vdW-DF [122], the ``opt" functionals (optPBE-vdW, optB88-vdW, and optB86b-vdW) where the exchange functionals were optimised for the correlation part [124], and the vdW-DF2 of Langreth and Lundqvist groups [125].

This method is available since the 5.2.12.26May2011 version of VASP for the calculation of total energies and forces. The stress calculation for the cell optimisation (ISIF=3) is available since the VASP 5.2.12.11Nov2011 version for spin unpolarised systems and VASP 5.3.1 for spin polarised systems.

N.B.: This feature has been implemented by J. Klimeš. If you make use of the vdW-DF functionals presented in this section, we ask that you cite the following paper:


J. Klimeš, D. R. Bowler, and A. Michaelides, Phys. Rev. B 83, 195131 (2011).


Correlation functionals

The method is invoked by setting

LUSE_VDW = .TRUE.

Moreover, the PBE correlation correction needs to be removed since only LDA correlation is used in the functionals. This is done by setting

AGGAC = 0.0000

The two tags above need to be used for all of the following functionals.


Exchange functionals

1) The GGA tag is further used to choose the appropriate exchange functional. The original vdW-DF of Dion et al uses revPBE, therefore the vdW-DF can be chosen by setting


GGA = RE

LUSE_VDW = .TRUE.

AGGAC = 0.0000


2) More accurate exchange functionals for the vdW correlation functional have been proposed in [124] and [126]. To use these functionals set:


GGA = OR

LUSE_VDW = .TRUE.

AGGAC = 0.0000


3) for optPBE-vdW,


GGA = BO

PARAM1 = 0.1833333333

PARAM2 = 0.2200000000

LUSE_VDW = .TRUE.

AGGAC = 0.0000




4)for the optB88-vdW functional, or


GGA = MK

PARAM1 = 0.1234

PARAM2 = 1.0000

LUSE_VDW = .TRUE.

AGGAC = 0.0000


5) for the optB86b-vdW functional.

In the vdW-DF2 functional the rPW86 exchange functional is used

GGA = ML

moreover, the vdW functional needs to be changed to the vdW2 correlation which requires only a change of a parameter:

Zab_vdW = -1.8867

Therefore to use vdW-DF2, set:


GGA = ML

LUSE_VDW = .TRUE.

Zab_vdW = -1.8867

AGGAC = 0.0000


An overview of the performance of the different approaches can be found for example in [124,125] for gas phase clusters and in [126] for solids.


126J. Klimeš, D. R. Bowler, and A. Michaelides, Phys. Rev. B 83, 195131 (2011).



Important remarks:

  • The method needs a precalculated kernel which is distributed via the VASP download portal (VASP -> src -> vdw_kernel.bindat) and on the ftp server (vasp5/src/vdw_kernel.bindat). If VASP does not find this file, the kernel will be calculated.

        This, however, is rather demanding calculation. The kernel needs to be either copied to the VASP run directory for each calculation or can be stored in a central location and read from there. The location needs to be set in routine PHI_GENERATE. This does not work on some clusters and the kernel needs to be copied into the run directory in such cases. The distributed file uses little endian convention and won't be read on big endian machines. The big endian version of the file is available from the VASP team.

  • There are no special POTCARs for the vdW-DF functionals and the PBE or LDA POTCARs can be used. Currently the evaluation of the vdW energy term is not done fully within the PAW method but the sum of the pseudo-valence density and partial core density is used. This approximation works rather well, as is discussed in [126], and the accuracy generally increases when the number of valence electrons is increased or when harder PAW datasets are used. For example, for adsorption it is recommended to compare the adsorption energy obtained with standard PAW datasets and more-electron POTCARs for both PBE calculation and vdW-DF calculation to assess the quality of the results.

  • The spin polarised calculations are possible, but strictly speaking the non-local vdW correlation is not defined for spin-polarised systems. For spin-polarised calculation the non-local vdW correlation energy is evaluated on the sum of the spin-up and spin-down densities.

  • The evaluation of the vdW energy requires some additional time. Most of it is spent on performing FFTs to evaluate the energy and potential. Thus the additional time is determined by the number of FFT grid points in the calculation, basically size of the simulation cell. It is almost independent on the number of the atoms in the cell. Thus the relative cost of the vdW-DF method depends on the ``filling" of the cell and increases with the amount of vacuum in the cell. The relative increase is high for isolated molecules in large cells, but small for solids in smaller cells with many k-points.

  • This feature has been implemented by J. Klimeš. If you make use of the vdW-DF functionals presented in this section, we ask that you cite the following paper: J. Klimeš, D. R. Bowler, and A. Michaelides, Phys. Rev. B 83, 195131 (2011).

 

3.VASP 手册解读: LVDW= .TRUE.     DFT-D2 method of Grimme

http://cms.mpi.univie.ac.at/vasp/vasp/DFT_D2_method_Grimme.html


DFT-D2 method of Grimme

LVDW= .TRUE. | .FALSE. (Available as of VASP.5.2.11)

Default: LVDW=.FALSE.


    Popular density functionals are unable to describe correctly van der Waals interactions resulting from dynamical correlations between fluctuating charge distributions.

   A pragmatic 【实际的, 注重实效的】method to work around this problem has been given by the DFT-D approach [120], which consists in adding a semi-empirical dispersion potential(Edisp) to the conventional Kohn-Sham DFT energy:


$.displaystyle E_{DFT-D} = E_{KS-DFT} + E_{.rm disp}.$(6.88)



In the DFT-D2 method of Grimme [121], the van der Waals interactions are described via a simple pair-wise force field, which is optimized for several popular DFT functionals. The dispersion energy for periodic systems is defined as


$.displaystyle E_{.rm disp} = -.frac{s_6}{2} .sum_{i=1}^{N_{at}} .sum_{j=1}^{N_{...
...,0}-{.bf r}^{j,{.bf L}}.vert^6} f(.vert{.bf r}^{i,0}-{.bf r}^{j,{.bf L}}.vert),$(6.89)



where the summations are over all atoms $ N_{at}$ and all translations of the unit cell $ {L}=(l_1,l_2,l_3)$, the prime indicates that $ i.not=j$ for $ {L}=0$, $ s_6$ is a global scaling factor, $ C_6^{ij}$ denotes the dispersion coefficient for the atom pair $ ij$, $ {r}^{i,{L}}$ is a position vector of atom $ i$ after performing $ {L}$ translations of the unit cell along lattice vectors. In practice, terms corresponding to interactions over distances longer than a certain suitably chosen cutoff radius contribute only negligibly to $ E_{.rm disp}$ and can be ignored. The term $ f(r^{ij})$ is a damping【阻尼;衰减】 function


$.displaystyle f(r^{ij}) = .frac{1}{1+e^{-d(r^{ij}/R_0^{ij}-1)}},$(6.90)



whose role is to scale the force field such as to minimize contributions from interactions within typical bonding distances. Combination rules for dispersion coefficients $ C_6^{ij}$ and vdW radii $ R^{ij}$ are


$.displaystyle C_6^{ij} = .sqrt{C_6^i C_6^j},$(6.91)



and


$.displaystyle R_0^{ij} = R_0^i+ R_0^j,$(6.92)



respectively. The global scaling parameter $ s_6$ has been optimized for several different DFT functionals such as PBE ($ s_6=0.75$), BLYP ($ s_6=1.2$), and B3LYP ($ s_6=1.05$).

The DFT-D2 method can be activated by setting LVDW=.TRUE. Optionally, the forcefield parameters can be controlled using the following flags (the default values are listed):


VDW_RADIUS= 30.0cutoff radius (Å) for pair interactions
VDW_SCALING= 0.75global scaling factor $ s_6$
VDW_D= 20.0damping parameter $ d$
VDW_C6= $ C_6$,...$ C_6$ parameters ($ Jnm^6mol^{-1}$) for each species defined in POSCAR
VDW_R0= $ R_0$,...$ R_0$ parameters (Å) for each species defined in POSCAR



The default values for VDW_C6 and VDW_R0 are compiled in Tab. 2. As the potential energy, interatomic forces as well as stress tensor are corrected by adding contribution from the forcefield, simulations such as the atomic and lattice relaxations, molecular dynamics, and vibrational analysis can be performed.

   The number of atomic pairs contributing to $ E_{.rm disp}$ and the estimated vdW energy are written in OUTCAR (check lines following the expression 'Grimme's potential'). The forces and stresses written in OUTCAR contain the vdW correction but the corrected energy should be read from OSZICAR (energies in OUTCAR do not contain the vdW term).

IMPORTANT NOTE: The defaults for VDW_C6 and VDW_R0 are defined only for the first five rows of periodic table of elements (see Tab. 2) - if the system contains other elements the user must provide the corresponding parameters.



Table 2: Parameters used in the empirical force-field of Grimme [121].
ElementC$ _6$R$ _0$ElementC$ _6$R$ _0$
Jnm$ ^6$mol$ -1$ÅJnm$ ^6$mol$ -1$Å
H0.141.001K10.801.485
He0.081.012Ca10.801.474
Li1.610.825Sc-Zn10.801.562
Be1.611.408Ga16.991.650
B3.131.485Ge17.101.727
C1.751.452As16.371.760
N1.231.397Se12.641.771
O0.701.342Br12.471.749
F0.751.287Kr12.011.727
Ne0.631.243Rb24.671.628
Na5.711.144Sr24.671.606
Mg5.711.364Y-Cd24.671.639
Al10.791.716In37.321.672
Si9.231.716Sn38.711.804
P7.841.705Sb38.441.881
S5.571.683Te31.741.892
Cl5.071.639I31.501.892
Ar4.611.595Xe29.991.881


119  A. Kiejna, G. Kresse, J. Rogal, A.De Sarkar, K. Reuter, and M. Scheffler, ``Comparison of the full-potential and frozen-core approximation approaches to density-functional calculations of surfaces'',
Phys. Rev. B 73, 035404 (2006).


120  X. Wu, M. C. Vargas, S. Nayak, V. Lotrich, and G. Scoles, ``Towards extending the applicability of density functional theory to weakly bound systems'',
J. Chem. Phys. 115, 8748 (2001).

121 S. Grimme,

``Semiempirical GGA-type density functional constructed with a long-range dispersion correction'',
J. Comp. Chem. 27, 1787 (2006).


网络问答:


A::::我要用VASP计算H2分子在表面的吸附,其中需要考虑范德华力,在INCAR中改如何设置参数呢?

Q:::

用DFT-D2方法,INCAR加入LVDW = .true.
用vdw-DF方法,加入LUSE_VDW = .TRUE.
这两种方法均有一堆其他参数需要设置,具体请参考手册。


A:::

用DFT-D2方法,INCAR只要加入LVDW = .TRUE. 就行了吧,像下列参数取默认值就行吧?
VDW_RADIUS = 30.0
VDW_SCALING = 0.75
VDW_D = 20.0
VDW_C6 =
VDW_R0 =


更多参考大牛文章:http://emuch.net/html/201203/4205083.html

从He2等的量子化学图形解读范德华作用(续)
作者: zhou2009 (站内联系TA)    发布: 2012-03-06



也用Gaussian09计算,用Multiwfn 2.3(含Gsgrid)进行数据处理,用Gaussview、Sigmaplot作图。当距离数值不标出单位时是Å,电荷和能量数值不标出单位时是 a.u.。

......

,范德华力与卡西米尔力都是实验证实客观存在的,它们也会有相互作用,但它们却不是同一个层次和范畴的。
       首先、从作用的形式和距离看,卡西米尔力是两块互相平行的中性导体板之间的吸引力,其间的距离数量级在微米,而范德华力是原子(分子)之间的,距离数量级在几埃米,两者之间相差上千倍!在卡西米尔力作用的距离,根本不会形成范德华力。而在范德华力作用距离,也已经谈不上卡西米尔力,因为这时形成卡西米尔力的平行板早崩塌了。
其次、从作用能和作用力看,范德华作用能数量级在kJ∙mol^-1,而卡西米尔作用能数量级在J∙mol^-1,卡西米尔力数量级在10^-7牛顿。
第三、从基本表达式看,范德华力F∝1/r^7,r是原子间距。卡西米尔力F∝1/a^4,a为平行板间距。
第四、从势能随间距变化曲线看,范德华作用有一个能量的最低点的平衡距离,如果距离再加大,能量会缓慢上升;如果从这平衡距离再缩小,能量会急剧升高。而卡西米尔作用,只有随间距缩小体系合拢受力永远增大这样一种趋势。
总的说来,范德华力与卡西米尔力是不同层次、不同性质的两种力。
       如果说范德华作用所形成的结合能比正常价键要小1~2个数量级的话,卡西米尔作用体系受力合拢能比范德华作用还要小2个数量级以上。
其实这二者在构成和性质上是完全不同的,本是没有可比性的。卡西米尔作用使体系受力合拢,只是数量上的变化,受外力合拢的体系本身材质不会有任何质的变化,这纯粹是个量子物理过程。而范德华作用所形成的结合由原子形成了分子、分子之间形成超分子,原子、分子发生了化学性质的变化,这种性质上的变化产生了范德华作用,这是个量子化学过程!
五、范德华作用量子化学图景的意义:
1、对范德华力的“电子的瞬间运动出现的瞬间偶极矩之间的作用” 的模型猜想,是仅用静电作用来解释范德华力,卡西米尔力更是在一定的格局下形成一种外力将材料挤拢,来解释范德华力,都脱离了量子化学。
这其实阻碍了对范德华力模拟、计算和研究。试想瞬间偶极矩如何度量、计算?并与实验测得的范德华作用距离、作用能相吻合?
        这使得我们通常对化学课题进行量子化学精确计算研究、建立了模型图像之后,一涉及到弱相互作用的范畴,就只是用范德华力一言而敝之,心里想着它是瞬间偶极矩作用,并常常就此止步了,对范德华作用不能通过计算形成一个具体的形成机理、图景。
          在量子化学看来,形成范德华力的双方电子是作为波在起作用的,而起作用的波是否合理,是用由它得到的电子密度分布安排是否合理来衡量的,而电子密度分布安排是要使得体系能量均衡(电负性均衡原理)、体系能量最低
      现在我们也动辄说电子密度、说电荷分布,但它们不是一种电荷模型设想,而是由量子化学计算的精细结果,电子是一个空间分布。范德华力是双方电子作为波相互作用之后的结果,包含了静电作用和电子相关,是内部起作用发生化学变化的结果,不是象卡西米尔力那样是在一定的格局之下,形成内外受力不均衡,由外力挤拢来形成结合力。也不是纯粹的静电引力作用,电子波的相互作用包含静电作用和电子相关,复杂得多。
2、现在我们看到了范德华作用的量子化学图景,它从成对电子的Pauli互斥开始,产生电子回避,形成了电子共享。我们不仅在惰性双原子准分子看到这种图景,而且对二体苯、二体叶绿素、生命分子的卷曲等等,都可以看到这种图景:电子排斥-回避-共享。
      其实精确的量子化学计算早已经完全包含了这个范德华作用的图景,只是我们因为我们囿于传统经典的“瞬间偶极矩”作用模型,没有通过电子密度差来具体观察、分析而已。
       可见只要选择的方法够好、基组合用,完全可以用来展现范德华作用的确切图景,研究范德华作用的具体过程和各种量值。吸附,催化,乃至化学反应等都是基于范德华作用而缘起、发生的。
3、现在有了以耦合簇为代表的量子化学方法,精确地考虑、计算了电子相关,从而也为我们树立了一个标杆。事实也正是这样,当Cybulski用CCSD(T)/aug-cc-pV5z方法计算了He2,Ne2,Ar2,He-Ne,He-Ar,Ne-Ar的基态势能曲线,预测了平衡键长、解离能等之后,许多新的方法、基组的建立都与它相对照,判断是否可用。
       现在对我们来说,对于一种量子化学方法,对于一组基组,以及它们是否配套,能不能用于研究我们的包含范德华作用的课题,我们可以在Gaussian09中用该方法/基组算一下惰性双原子准分子,并作出电子密度差,来与CCSD/aug-cc-pVTz算的密度差对照一下,就可以看出我们采用的方法/基组是否合用,合用的程度。并分析不合用的原因。
比如用CCSD/aug-cc-pVTz与CCD/aug-cc-pVTz作的电子密度差基本没有区别。但是对于MP4SDQ/cc-pVTz与MP4SDQ/aug-cc-pVTz,它们完全不相同,只有后者(加了aug)才能与CCSD/aug-cc-pVTz基本相符。
4、        要使用电子密度差方法来查看化学图景,这使得我们似乎看见了化学。
         现在用它研究各种化学体系的范德华作用。主要观察什么部位是排斥的中心,由于排斥电子又回避到什么部位了,在什么部位形成了什么样的共享电子图景。当然,如果涉及的是催化、化学反应,除了分析体系的范德华作用本身外,还要对体系优化结果,考察范德华作用引起的其他正常键的松动、变化,这里正是催化、化学反应的本源。
5、当然,根据耦合簇为代表的多种量子化学方法,精确地考虑、计算电子相关,得到的电子密度差图景,从而跟踪到的范德华作力形成的过程及本质,还需要进一步有更多的数据、事实来深入探讨和论证。特别是所谓共享电子,它有多少?是否足以维系着体系形成范德华作用的结合能,需要进一步数据确证。
6、物理有机,乃至化学,历来有多少经验的猜想、假说、原理,现在需要我们用严格的量子化学计算、图像来重新考察一下啊! (2012/3/5)
 




https://blog.sciencenet.cn/blog-567091-731874.html

上一篇:Phonopy test:反铁磁、铁磁、虚频、k-mesh密度
下一篇:为了忘却的纪念:Fireworks使用简介及其他
收藏 IP: 128.84.126.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-9-15 19:14

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部