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

博文

结构优化效率与计算参数设置

已有 26575 次阅读 2013-6-1 23:14 |个人分类:结构优化|系统分类:科研笔记| 计算, 结构

0.加快收敛速度之EDIFF参数设置


EDIFF-tag EDIFF= [real]

Default :

EDIFF=$ 10^{-4}$

  Specifies the global break condition for the electronic SC-loop. The relaxation of the electronic degrees of freedom will be stopped if the total (free) energy change and the band structure energy change ('change of eigenvalues') between two steps are both smaller than EDIFF.

    For EDIFF=0, NELMelectronic SC-steps will always be performed.

Mind: In most cases the convergence speed is exponential. So if you want the total energy significant to 4 figures set EDIFF=$ 10^{-4}$. There is no real reason to use a much smaller number【在总能当中电子部分的能量起决定作用,总能一定意义上即电子部分的能量,只有计算声子谱才考虑到了晶格振动部分的能量】



EDIFFG-tag

EDIFFG= [real]


Default:

EDIFFG=EDIFF*10

   EDIFFG defines the break condition for the ionic relaxation loop. If the change in the total (free) energy is smaller than EDIFFG between two ionic steps relaxation will be stopped.

    If EDIFFG is negative it has a different meaning: In this case the relaxation will stop if all forces are smaller than $ .vert$EDIFFG$ .vert$. This is usually a more convenient setting.

  EDIFFG might be 0; in this case the ionic relaxation is stopped after NSW steps.

EDIFFG does not apply to MD-simulations.

网上解读:

    EDIFF=1E-5是电子自洽步(SCF)的收敛标准,表示0.01meV/cell;
    EDIFFG表示离子弛豫的收敛标准,若为正值,表示以能量为收敛标准,意义同上,默认值为EDIFF的10倍;若为负值,表示以离子受力作为收敛标准,单位eV/A,这个标准不能太高,否则很难收敛。
    一般情况下,你可以先设置EDIFFG为正值,跑完后用p4vasp查看一下,能看到能量和力随弛豫步的变化(跟MS类似),能量最低的时候力也基本上就是最小的了,而且对稍微复杂一点的体系,力很难收敛到0.001eV/A。

//////////////////////

对于大的体系,还是分步来优化好,先把EDIFF EDIFFG弄大一点 ,设置POTIM =0.1 那个NSW的步数也设小一点,然后把EDIFF EDIFFG弄小一点 设置POTIM =0.2 NSW=1000 把那个IALGO去掉,进行第二次优化

/////////////

算法也决定你计算的准群度呀,我一般用48做最开始优化,后来就用默认的,k点用几个你自己可以测试一下,不能只看晶胞大小

//////////



1.截断能对计算速度的影响有多大?

  A:书上只是说,截断能必须进行熟练测试。所以我们就想当然的对截断能进行收敛测试。但是,为什么不直接选择一个很大的截断能呢?这样一定能保证收敛。而且还省去了收敛测试的步骤。也没有表明截断能太大会影响计算的速度。

     Q: 谁说截断能不影响计算速度,截断能取得大,电子波函数的平面波展开中就要取更多的项。平面波越多,计算量就越大。

  截断能是表示平面波展开后取到多大能量的平面波,对于高能部分,展开后所占的比例非常小,而且影响计算速度,所以并不是截断能越大越好。这是一个两相竞争的问题。一般用POTCAR中的ENMAX即可,若做静态等计算,可适当调大些。

Q2:

控制参与计算的平面波的个数,截断能约大,所参与的平面波就越多,计算的精度当然会随之提高,然而,在超过一定的数值后,计算的结果没有实质的变化,反而计算的时间却增加了

Q3:这个一定要测试的。还是要把计算理论好好看看的。

Q4:

一般,认为截断能越大,k点设置越多,收敛精度越高,收敛会越慢,但是准确度不一定越高。因此有必要做一些此类参数的测试。一般会这些参数大概会有个能量收敛的趋势。就是设置到某个临界参数,能量几乎变化很小。而参数再增加,只是增加计算量,而准确度没有变化。
对于比较大很难优化的周期体系,一般建议是采用 1 1 1的k点先做一个优化,取得一个合理的初步结果,然后采用大一点的k点设置进行计算。可能结构优化速度会快很多。一个朋友计算过一个大体系,发现结合两者的方法(1周)比直接采用后者的方法(4周)节省了3周。
另外也可以采用先用一个低等的精度做一个预优化,然后高点精度的参数设置进行优化。

做结构优化是做计算的必做任务。大家有什么好的经验一起讨论下,一起进步下:-)

Q6:

优化时可以用ismear=0,最后静态的时候再用-5,这样有些错误就不会发生。
楼主说的分级优化是绝对可行的,vasp说明书也说了优化完后要在优化一次,然后再做static计算。
实际上应该优化到一个很高的精度再做静态计算会比较合理。

Q7:

关于截断动能的优化,为什么要设置不同的ENCUT进行计算,以得到一个所谓较好的值(判断依据好像是总能变化在0.001eV)。是单纯从节省计算时间角度考虑么?否则的话直接取相对较大的ENCUT(ENMAX的1.3倍)不是可以省去优化ENCUT这一步了吗?但是假如为了优化ENCUT做了6次计算,而且有的ENCUT设置到很大,那么为了优化ENCUT所花费的时间好像和所说的以适当ENCUT值做后面的性质计算所节省的时间也差不多吧?是不是后面所做的结构优化或者性质计算里参数改变的计算次数很多?或者有什么更细节的步骤,希望大家能够不吝赐教。

Q8:

节省计算时间。
后面需要用到的次数可能还会很多,科研不是做一天两天的事情。完全有必要找一个合理的截断能。否则写文章的时候,你直接写个很大的截断能(准确性没问题),但是会让人怀疑你对计算的理解,其他的结果都可能被怀疑了。

Q9:

对于具有周期性的固体而言,由于有布洛赫原理,说明了周期波函数描述的正确性,从而使平面波的描述是一个比较合理和最为理想的方案。截断能是指k空间的截断,作为基组,原则上需要完备基组,这样要求的是截断能尽可能地大,以至于使基组近乎完备,但是过高的截断能对结果的改善并不太明显,所以作为一般的计算,采取一定的截断能可以使计算的精度和计算的可行性有一个比较合理的平衡。所以选用平面波基组对于晶体的计算时一个好的选择,采取合适的截断能可以计算出比较接近真实情况的性质。


A2: The Choice of a Basis Set - Plane Waves


The Kohn-Sham orbitals may be represented in terms of any complete basis set. Many choices are possible including atomic orbitals, Gaussians, LAPW  and plane waves, the basis set we use in practice. The use of a plane wave (PW) basis set offers a number of advantages, including the simplicity of the basis functions, which make no preconceptions regarding the form of the solution, the absence of basis set superposition error, and the ability to efficiently calculate the forces on atoms  

In general, the representation of an arbitrary orbital in terms of a PW basis set would require a continuous, and hence infinite, basis set. However, the imposition of periodic boundary conditions allows the use of Bloch's Theorem  may be written



where the sum is over reciprocal lattice vectors   is a symmetry label which lies within the first Brillouin zone. Thus, the basis set for a given   will be discrete, although in principle it will still be infinite. In practice, the set of plane waves is restricted to a sphere in reciprocal space most conveniently represented in terms of a cut-off energy,  such that for all values of  used in the expansion




   Thus, the convergence of the calculation with respect to basis set may be ensured by variation of a single parameter, . This is a significant advantage over many other basis set choices, with which calculated properties often show extreme sensitivity to small changes in basis set and no systematic scheme for convergence is available.

The choice of periodic boundary conditions is natural in the case of bulk solids which exhibit perfect translational symmetry.

Q1:

根据前面的讨论知道,在程序的实际计算中,采用平面波基组时,每个本征值对应的波函数都是选用的平面波的线性叠加。平面波的个数N取决于截断能的大小Ecut,具体数值可以根据9楼公式2.13得到。计算时所要做的是优化每个平面波的系数,迭代求解,直至自洽,这样得到了每个本征值所对应的波函数中N个平面波叠加的一套系数和相应的本征值。(题目一:假设一个单胞的晶格常数分别是a=b=c=10 埃,截断能是400eV,那么每个波函数需要多少个平面波叠加得到?

Q2:k电荷G没有直接关系,K点是体现周期性的,根据玻恩-冯卡门条件得到;G是根据截断能得到的,一旦截断能确定G是确定的。之所以式子中出现K+G,目的是对每个K点的波函数的展开都需要G个平面波的展开。


Q3:

在式子2.13中
0.5(k+G)^^2  <= ENCUT
一个G对应于一个平面波,那是不是说k越大,G点的数目也就小了,也就是说对应的平面波数,也少了? 下面的结果是我得到的平面波数与vasp中OUTCAR中的不一致,不知道怎么回事,当然很有可能是程序的问题。


A3:从量子化学到固体理论-内在的过渡

    不明白K点到底是什么?又为什么波函数要在K点展开?“波函数在某点展开”又是什么意思?请版主明示。

    按照我的理解,
分子中电子的波函数等于电子在孤立原子中波函数的线性组合,此即分子轨道理论而孤立原子中的电子的波函数即为基组

     为什么选用孤立原子中的电子的波函数作为基组进行线性组合呢?因为电子的量子理论基础,是说我们并不知道电子处于哪个态(比如1s, 2s, 2p, 3s, 3p, 3d,...态),我们只知道它是它所有可能态的线性组合,因此需要找到所有的可能态进行组合,这就是所谓的基组要“完备”,即找到所有的可能的态    这里隐含的意思就是说基组所代表的就是电子的一个“态”。

     其次,基函数还要“正交”,即基函数所代表的态之间彼此独立,因为函数正交在数学上是说两函数乘积的积分为0,物理上讲是说任意两态之间交叠的概率为0. 因此,我们可以选用一套表示电子在各种“态”的“正交“,”完备”的基函数作为起点,通过求解薛定谔方程,来优化任何多电子体系中每个电子的波函数,最后总波函数就是所有电子波函数的乘积。这就是整个量子化学的起点,HF理论。

     这里我们试图理解的是,为什么选用平面波来做第一原理计算,但是中间还隔着个密度泛函理论。我觉得简单说,有人(H-K)证明了电子体系基态的各种动能,作用能与密度之间存在一一对应关系,那么可以以密度为变量,通过优化电子概率密度来找到体系基态能量。

     像电子动能,电子-核吸引,电子电子排斥,这些项的能量与密度之间的关系比较简单,库仑作用嘛,遗憾的是不知道这个所谓的交换相关能与密度之间的物理和数学上的“一一对应”是怎么个对应法,也就是说不知道其泛函,因此提出了很多泛函类别。

     但是这里要讨论的问题并不在于这里,而是平面波基组为什么能表示周期体系的电子量子”态“,又如何表示。首先密度泛函理论说,因为我只关注密度到能量的表示是否正确,因此我不care密度是怎么算来的,原则上讲,任何完备,正交的基函数,比如各种多项式,都可以用来计算电子的密度。但这只是原则,实际上非周期的分子体系中,电子的波函数还是通过孤立原子的”态“进行线性组合以计算密度的,即DFT在密度这一块和HF理论完全一致,没啥差别,只是在从密度到能量的计算上不一致。
      基于上面的理解,我觉得
平面波被选为周期体系的电子”态“的数学表示,主要的优势是其具有周期性,这就是Bloch theorem吧,就是说Bloch的功劳,是找到了一套数学函数,具有奇妙的周期性,完备性和正交性,可以用来表示无限平移周期体系中的电子。

     但是,从这里,一直到波函数在倒空间中的K点上展开,我就始终参不透了。所有讲述这方面理论的材料就像讲天书一样不说人话【经典! 人话就是用自己的话来诠释!】。恳请版主指点。特别是以下两个疑问:

”K点到底是什么?又为什么波函数要在K点展开?“波函数在某点展开”又是什么意思?“

       相信K点在这里也代表电子的各个”态“,比如1s, 2s, 2p, 3s, 3p, 3d,...至少对应这样的物理意义,那为什么要叫什么”波矢“呢?波函数在K点展开,应该是在这里进行线性组合;HF理论通过用不同的Slater指数函数参数来表示离核远近的电子概率分布不同,以形成s, p, d, f那些奇异的轨道形状,平面波函数貌似都是那一个形式,它是如何做到的呢?


2.  关于结构弛豫设置及收敛判据选择的一点体会

摘录自:http://emuch.net/bbs/viewthread.php?tid=2512497

(1)收敛判据的选择

      结构弛豫的判据一般有两种选择:能量和力。这两者是相关的,
理想情况下,能量收敛到基态,力也应该是收敛到平衡态的。但是数值计算过程上的差异导致以二者为判据的收敛速度差异很大,力收敛速度绝大部分情况下都慢于能量收敛速度。这是因为力的计算是在能量的基础上进行的,能量对坐标的一阶导数得到力。计算量的增大和误差的传递导致力收敛慢。

   到底是以能量为收敛判据,还是以力为收敛判据呢?关心能量的人,觉得以能量为判据就够了;关心力相关量的人,没有选择,只能用力作为收敛标准。对于超胞体系的结构优化,文献大部分采用Gamma点做单点优化。这个时候即使采用力为判据(EDIFFG=-0.02),在做静态自洽计算能量的时候,会发现,原本已经收敛得好好的力在不少敏感位置还是超过了结构优化时设置的标准。这个时候,是不是该怀疑对超胞仅做Gamma点结构优化的合理性呢?是不是要提高K点密度再做结构优化呢?

   在我看来,这取决于所研究的问题的复杂程度。我们的计算从原胞开始,到超胞,到掺杂结构,到吸附结构,到反应和解离。每一步都在增加复杂程度。结构优化终点与初始结构是有关的,如果遇到对初始结构敏感的优化,那就头疼了。而且,还要注意到,催化反应不仅与原子本身及其化学环境有关,还会与几何构型有关。气固催化反应过程是电子的传递过程,也是分子拆分与重新组合的过程。如果优化终点的构型不同,可能会导致化学反应的途径上的差异。仅从这一点来看,第一性原理计算的复杂性,结果上的合理性判断都不是手册上写的那么简单。

对于涉及构型敏感性的结构优化过程,我觉得,以力作为收敛判据更合适
而且需要在Gamma点优化的基础上再提高K点密度继续优化,直到静态自洽计算时力也是达到收敛标准的。

(2)结构优化参数设置

       结构优化,或者叫弛豫,是后续计算的基础。其收敛性受两个主要因素影响:初始结构的合理性和弛豫参数的设置。

初始结构

    初始结构包括原子堆积方式,和自旋、磁性、电荷、偶极等具有明确物理意义的模型相关参数。比如掺杂,表面吸附,空位等结构,初始原子的距离,角度等的设置需要有一定的经验积累。DFT计算短程强相互作用(相对于范德华力),如果初始距离设置过远(如超过4埃),则明显导致收敛很慢甚至得到不合理的结果。

比较好的设置方法可以参照键长。比如CO在O顶位的吸附,可以参照CO2中C-O键长来设置(如增长20%)。也可以参照文献。记住一些常见键长,典型晶体中原子间距离等参数,有助于提高初始结构设置的合理性。实在不行,可以先在小体系上测试,然后再放到大体系中算。

弛豫参数

弛豫参数对收敛速度影响很大,这一点在计算工作没有全部铺开时可能不会觉察到有什么不妥,反正就给NSW设置个“无穷大”的数,最后总会有结果的。但是,时间是宝贵的,恰当的设置3小时就收敛的结果,不恰当的设置可能要一个白天加一个黑夜。如果你赶文章或者赶着毕业,你就知道这意味这什么。

结构优化分电子迭代和离子弛豫两个嵌套的过程。

电子迭代自洽的速度,有四个响很大的因素:初始结构的合理性,k点密度,是否考虑自旋和高斯展宽(SIGMA);

离子弛豫的收敛速度,有三个很大的影响因素:弛豫方法(IBRION),步长(POTIM)和收敛判据(EDIFFG).

一般来说,针对理论催化的计算,初始结构都是不太合理的。因此一开始采用很粗糙的优化(EDIFF=0.001,EDIFFG=-0.2),很低的K点密度(Gamma),不考虑自旋就可以了,这样NSW<60的设置就比较好。其它参数可以默认。

经过第一轮优化,就可以进入下一步细致的优化了。就我的经验,EDIFF=1E-4,EDIFFG=-0.05,不考虑自旋,IBRION=2,其它默认,NSW=100;跑完后可以设置IBRION = 1,减小OPTIM(默认为0.5,可以设置0.2)继续优化

优化的时候让它自己闷头跑是不对的,经常看看中间过程,根据情况调节优化参数是可以很好的提高优化速度。这个时候,提交两个以上的任务排队是好的方式,一个在调整的时候,下一个可以接着运行,不会因为停下当前任务导致机器空闲。

无论结构优化还是静态自洽,电子步的收敛也常常让新手头痛。如果电子步不能在40步内收敛,要么是参数设置的问题,要么是初始模型太糟糕(糟糕的不是一点点)。【如何查看电子步的收敛!!!】

     静态自洽过程电子步不收敛一般是参数设置有问题。这个时候,改变迭代算法(ALGO),提高高斯展宽(SIGMA增加),设置自洽延迟(NELMDL)都是不错的方法。对于大体系比较难收敛的话,可以先调节AMIN,BMIX跑十多步,得到电荷密度和波函数,再重新计算。实在没办法了,可以先放任它跑40步,没有收敛的迹象的话,停下来,得到电荷密度和波函数后重新计算。一般都能在40步内收敛。

对于离子弛豫过程,不调节关系也不大。开始两个离子步可能要跑满60步(默认的),后面就会越来越快了。

总的说来,一般入门者,多看手册,多想多理解,多上机实践总结,比较容易提高到一个熟练操作工的水平。

   
如果要想做到“精确打击”,做到能在问题始发的时候就立刻采取有效措施来解决,就需要回归基础理论和计算方法上来了。

高手啊高手,到底是什么样的境界?

(3)优化结果对初始结构和“优化路径”的依赖

       原子吸附问题不大,但是小分子吸附,存在初始构型上的差异。slab上水平放置,还是垂直放置,可能导致收敛结果上的差异。根据H-K理论,理想情况下,优化得到的应该是全局最小,但在数值计算的时候可能经常碰到不是全局最小的情况。实际操作中发现,多个不同初始结构优化收敛后在能量和结构上存在一定差异。

为了加快收敛速度,特别是对于表面-分子吸附结构,初始放松约束,比如EDIFF=1E-3,EDIFFG=-0.3,NSW=30可能是很好的设置。但是下面的情况应当慎重:
EDIFF=1E-3;
EDIFFG=-0.1;!或者更小
NSW=500;!或者更大

电子步收敛约束较小,而离子步约束偏大,离子步数又很多,这种情况下,可能导致的结果是结构弛豫到严重未知的区间。

再在这个基础上提高约束来优化,可能就是徒劳的了——结果不可逆转的偏向不正常的区间。

好的做法,是对初始结构做比较松弛的约束,弛豫离子步NSW应该限制在一个较小的数值内。EDIFF=1E-3的话,EDIFFG也最好是偏大一些,如-0.3而不是-0.1. 这样可以在较少的步数内达到初步收敛。

对于远离基态的初始结构,一开始在非常松弛的约束下跑若干离子步,时间上带来的好处是很大的。对于100个原子的体系用vasp做Gamma点优化,如果一开始就是正常优化( EDIFF=1E-4,EDIFFG=-0.02)设置,开始十个离子步可能都要花上几个小时。如果这个时候才发现输入文件有错误,那下午的时间就白费了,顺便带上晚上机器空转。

所以,我习惯的做法,是在初始几步优化后,会用xcrysden 检查一下XDATCAR中的数据,用xdat2xyz.pl生成movie.xyz,然后看看弛豫过程是不是按照设想的那样。后续过程跑完一个收敛过程,就再检查一下movie.xyz。如此这般,才放心的展开后续计算。

结构优化是个烦心的事情,这是所有可以写入论文中的甜蜜内容的前戏。如果到后来处理数据画出漂亮的图形来才发现结构优化错了,那时就会感叹手册上提醒的:test,test,test,。。。。是多么的高明。


本文来自: 小木虫论坛 http://emuch.net/bbs/viewthread.php?tid=2512497&fpage=1&view=&highlight=&page=5

(4)目的导向的结构优化

结构优化到这个阶段,是高级的了。为了得到特定结构,或者为了验证某些猜想,需要设计合理的初始结构,然后在这个基础上小心优化,比如POTIM=0.1跑几步看看,然后修改优化参数。

我遇到过的一件跟结构优化关系很大的算例是CeO2氧空位结构电子局域的问题(http://emuch.net/bbs/viewthread.php?tid=2954558)。按照一般方式(从优化好的bulk建slab模型,然后优化)得到一个O空位留下的两个电子均匀局域到O次外层三个Ce原子上,得到空位形成能2.34eV.经高人指点后,调节空位附近O原子位置,打破对称性后重新优化,两个电子完美的局域到两个Ce原子上了。并且空位形成能降低到2.0X eV。从这个例子可以看到,结构优化存在不少技巧的,这些技巧建立在研究者对模拟对象的物理意义的理解上。对物理图像的直观深入理解,才能做好模型预设,在此引导下才可能有目的的优化出不比寻常的结果。

目前第一性原理理论中的交换关联泛函部分包含经验参数。考虑这一点对优化结果的影响也很有意思。比如有专家提到,DFT+U参数对某些结构的收敛终态构型有影响。构型的变化可能影响表面反应过程。基于这一点,一个好的计算研究可能就出来了。

真实过程总是复杂多变的。无论何种模拟,估计都可以找到一些试验现象来验证。但是到底应该如何评判模拟结果,如何从第一性原理研究中得出有意义的结论需要很好的洞察力。这样的模拟不见得就必须建立的试验的基础上,完全凭空设计的模型有可能更能优美的解释本质。

今天是4月1日,是我注册小木虫3周年纪念日,上午领了60个金币大礼包。祝大家万事如意,学业有成。


Q1:

初略计算(示意设置,请勿复制):

ISTART =0; ICHARG =2;EDIFF=1E-3;EDIFFG=-0.2,ISPIN=1;NSW=40

再计算(示意设置,请勿复制):
cp CONTCAR POSCAR
ISTART =1;ICHARG=1;EDIFF=1E-4;EDIFFG=-0.02;ISPIN=2;NSW=200

Q2:

G点是对称性最高的点,简并度最高(这个,我不是很确定呢,呵呵),因而取一个K点的时候,就取Gamma点最具代表性。而且取Gamma点,程序计算速度最快。

对于分子,团簇等本身就是非周期性结构,一般都是推荐用单K点优化。这源于其物理意义。但是在用周期性结构来模拟分子和团簇的时候,我不敢说不用加密K点。

Q3:

>xdat2xyz XDATCAR

!上面这个命令生成movie.xyz

>xcrysden --xyz movie.xyz

上面这个命令图形化打开movie.xyz

这个是vtstool code中的命令。去这里下载下来:http://theory.cm.utexas.edu/vtsttools/


Q3:

1.POTIM在优化时是粒子移动的一个系数,也就是决定下一步移动的距离
2.ISTART程序会自己判定,即如果没有WAVECAR,但是令ISTART=1也是不起作用的
3.EDIFF是能量判据
EDIFFG设定为负值的话是力的判据


Q5:

截断能是在计算赝势中把价态电子和核(芯态)电子分开(cut off)的能量。赝势的目的就是得到一个势来最好地描述其对价态电子的行为, 因此如何选择价态与芯态对赝势行为有本质的影响.  在确定了价态与芯态之后, 我们需要选择核区半径(截断半径)来截断核区电子的波函数,也就是所谓match radii, 这个半径选取直接影响生成赝势的行为.

有网友提出疑问:
1第二行的“其”指什么?
2、“因此如何选择价态与芯态对赝势行为有本质的影响. 在确定了价态与芯态之后, 我们需要选择核区半径(截断半径)来截断核区电子的波函数,”不知道是先确定赝势还是价态电子
3、核区半径(截断半径)与核区电子的波函数有什么关关系?




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

上一篇:科研小记:跟踪自己关心的领域的方法(摘录)
下一篇:结构优化效率与计算参数设置II-VASP手册解读
收藏 IP: 67.255.15.*| 热度|

0

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

数据加载中...
扫一扫,分享此博文

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

GMT+8, 2024-5-19 08:37

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部