0%

Hall Petch关系是材料学中描述多晶材料最为重要的一个相关关系,它反映了材料强度与晶粒的尺寸的开方成反比。同时材料在极小晶粒尺寸下表现出的反常Hall Petch 关系也一度成为材料学家的研究热点。关于反常Hall Petch关系的研究发现屡现于高水平期刊,本文将不做更多介绍。从Hall Petch关系的第一次报导开始,已经有大量的文章针对各种材料体系(纯金属,金属间化合物,多相合金等等)进行了相关研究。然而,从Hall Petch 关系提出至今,依旧没有出现一个完美简洁的理论可以解释这个现象。研究人员为了解释这个现象已经提出了大量的理论,在接下来的篇幅中,作者将从中提取了几个典型的理论做简要介绍

1951年,Hall(文献1)发表的论文中首次报导了钢中存在的Hall Petch关系:

$\sigma=\sigma0 + k{HP}d^{-1/2} $

其中$\sigma$表示流变应力,d表示晶粒尺寸,k是Hall Petch 系数,其对不同材料,温度和应变条件都非常敏感。

与此同时,他也在文章中首次提出了试图解释Hall Petch关系的pile up(位错阻塞)模型。

该理论认为材料的强度来自于外力作用下导致的晶粒内部位错源开动,位错源启动后形成的位错林被晶粒之间的界面(晶界)阻碍后引起了材料的应力强化。简单起见,我们不妨假设两个近似六边形的晶粒以晶界AB接触,左边晶粒的中心O点存在一个位错源,在外加应力的作用下位错源沿着某一滑移面释放出一系列位错。随着第一个发射的位错冲撞到晶界AB上,整个位错列队被禁锢在位错源与晶界之间。由于晶界墙的阻力和新产生位错的弹性应变,整个位错列最后会形成一个稳定分布,假设晶粒的平均尺寸为d,则位错列中位错的数目为:

$n=\frac{(1-v)d\tau_e}{2\mu{b}}$

其中$\tau_e$表示位错的有效剪切应力,它可以表示为:

$\tau_e=\tau-\tau_i$

显然,它实际上是在加载应力中扣除了弹性应力部分。试想如果你要产生一根位错,首先晶体要发生塑性变形,因此塑性变形前的弹性部分并不是我们应该要考虑的部分。既然我们已经知道位错列中位错的数目,如果我们粗略地认为,所有的位错作用到晶界上的应力都是等效的,那么晶界上受到的应力$\tau_p$为

$\tau_p = n\tau_e$

Hall认为,材料屈服发生于位错列冲破晶界的时候,也就是说 $\tau_p$大于临界值$\tau_c$。如果我们取等号,再将前面的n代入,我们就可以得到:

$\tau = \tau_i +\surd{\frac{2\tau_c\mu{b}}{1-v}} \times d^{-1/2}$

该公式看起来非常合理简洁,-1/2次方的关系也出来了。但是如果材料的屈服依靠的是位错像裂纹一般冲破晶界,我们应该在实验中大量发现,但实际上,只有少数几篇文献报导过相关现象。更重要的是,依据该理论计算的屈服强度跟实际数据还是有较大出入。

在pile up模型的基础上,Cotrell(文献2)认为,材料的屈服并非来自于位错林冲破晶界,而是由于位错林在晶界的终端的应力场诱发了邻近晶粒中的位错源(如下图的R处)的开动。这样子需要的 临界应力要小于pile-up模型。

晶粒内部出现位错源固然可行,但是另一方面,晶界本身就是天生的位错源,在晶粒内部没有经受剧烈变形的情况下,晶界上的位错形核位置要远远多于晶粒内部。Li (文献3)于1963年提出了boundary source(晶界位错源)模型。为了方便起见,我们不妨假设晶粒是边长为d的立方体。假设晶粒的每一个面上的台阶密度为m,假设每一个台阶都有一定几率发射位错。进一步简化,我们不妨认为每个台阶都会发射位错,则一个面的总位错数目为$md^2$。每一个晶粒与毗邻晶粒公用一个晶面,则一个晶粒可以认为有三个晶界面属于这个晶粒。所以:

我们可以认为一个晶粒有3$md^2$个位错。这些位错在外加应力作用下继续向晶内扩散,形成位错林。位错林的密度为:
$\rho = 3md^2 / d^3 = 3m/d$
另一方面,晶体材料的位错密度和强度之间可以由著名的Taylor公式来描述:
$\sigma = \sigma_0 + M\alpha\mu{b}\surd\rho$
代入前者即得:

$\sigma = \sigma_0 + M\alpha\mu{b}3^{1/2}m^{1/2}d^{-1/2}$

我们再一次获得了Hall Petch关系。在原文中也有说明这个方式计算的理论屈服强度与pile up模型基本一致。目前介绍的这两个模型,前者认为位错完全起源于晶粒内部,后者则认为位错完全起源于晶界,两者分别是实际情况下的两个极端。但是殊途同归,两个极端最后得到了相同的结论和近似的屈服强度。不过,这两个模型的局限性也是明显的,实验中Hall Petch曲线的斜率对晶粒结构和化学成分非常敏感,但是上述模型中并未考虑这些实际因素。

值得注意的是,在2002年的一篇文章4中,该文作者成功在MD模拟中发现了晶界处位错的形核和发射过程。

(未完待续)

1 Hall, E. O. The deformation and ageing of mild steel: III discussion of results. Proceedings of the Physical Society. Section B 64.9 (1951): 747.
2 A.H. Cottrell,: The Mechanical Properties of Matter, Wiley, New York (1964)
3 Li, James CM. Petch relation and grain boundary sources. Transactions of the Metallurgical Society of AIME 227.1 (1963): 239.
4 Van Swygenhoven, H., P. M. Derlet, and A. Hasnaoui. Atomic mechanism for dislocation emission from nanosized grain boundaries. Physical Review B 66.2 (2002): 024101.
本文经ponychen授权发布,版权属于ponychen。

基于DFT+U研究应力调控CeO2(111)催化分解水制氢的活性

研究背景

氢能清洁无污染,是最具潜力的能量载体。使用太阳能、风能和水能转换得到的清洁电能电解水制氢是人来未来最具发展潜力的、可持续的、绿色环保的制氢技术。传统低温碱性电解水制氢技术的效率受限于液体电解质较低的传质过程和氧电极较慢的反应动力。相比于液体传质技术,新兴发展的高温固体燃料电池 (SOEC)技术以固体作为传质媒介,充分利用工业废热提供反应所需的热能,具有极高的传质效率和水电解效率。然而人们对SOEC制氢反应机理尚缺乏深刻认识。计算模拟技术是揭示能源材料领域中多相催化反应过程及机理的重要手段,可以辅助了解SOEC电解水的关键反应步骤,以开发、低廉且稳定的催化剂用于增强氢电极的稳定性以及提高制氢反应在燃料电池中的效率。兼具优异的电子和离子传质能力的氧化铈(CeO2)是SOEC水分解重要的催化剂和稳定剂。CeO2优越的催化活性归因于Ce3+与Ce4+的氧化还原循环。基于密度泛函理论,该工作使用Hubbard-UUeff=4.5 eV)修正的方法来描述Ce 4f轨道电子的强关联作用,研究表面应力对CeO2(111)催化制氢的影响。

研究内容

水在CeO2表面分解制氢包括三个基本反应步骤:(1)CeO2(111)表面生成一个氧空位;(2)H2O吸附于氧空位附近最终形成表面OH-;(3)表面OH-分解形成H2并留下一个晶格氧,如Figure 1所示水分解过程中形成的反应中间体以及最终生成H2的反应过程。该工作主要研究了三种可能的反应路径,即不同OH- 覆盖率下形成H2的过程如Figure 1所示。通过研究应力对三种不同反应路径中反应中间体形成的影响,以及应力对各反应步骤反应能垒的调控,来探究应力如何增强CeO2催化制氢的性能。

首先,研究发现拉伸应变有效增强了氧空位和OH-等中间体的稳定性,随之增加了各反应步骤的反应能垒。与之相反的是压缩应变可以通过减弱OH-与表面的结合能,从而降低生成H2的能垒。通过计算不同反应温度下中间体的自由能可以获得沿着2H,4H和5H的反应路径能量图。研究发现不管处于何种应力条件下生成氧空位和氧空位的扩散,以及OH-的形成不需要很高的能量势垒,而OH-分解形成H2是水分解制氢反应路径中能垒最高的一个步骤。然后,通过对比不同反应温度下,研究发现在温度低于1000 K时,沿着5H路径可以最快速地水分解生成H2 如图2所示,且最佳反应路径不受应力的影响。然而,当温度升高到1200 K时,应力有效地调控了制氢的反应路径,即拉伸应力下最优路径沿着5H,而在压缩应力下最优反应路径沿着2H或者4H。此外,研究还发现反应速度随着压缩应力的增大而增快,且当压缩应力大于3%时,反应速度获得显著提升。因此,计算模拟理论研究显示压缩应力(>3%)有利于增强CeO2(111)催化水分解制氢的反应活性。

图1. 2H,4H和5H覆盖的CeO2(111)表面生成H2的反应路径示意图。黄色, 红色,蓝色,白色,黑色和灰色原子球分别代表Ce,表面O,次表面O,H,次表面氧空位和表面氧空位。

图2不同反应温度与应力条件下最优反应路径生成H2的反应速度对比图。

作者简介

本文第一作者,武甜甜为丹麦科技大学能源转化与储存系博士后研究员。

Google Scholar:

https://scholar.google.com/citations?user=rm84jfsAAAAJ&hl=en

指导老师为:Tejs Vegge教授。

课题组主页:https://www.dtu.dk/english/service/phonebook/person?id=5334&tab=1。

本文经作者授权发布,版权属于第一作者。

理论模拟导向的催化剂理性高效设计策略用于醇的选择性胺化

第一作者:王涛

通讯作者:Philippe Sautet

研究背景

胺是工业上生产药物、农用化学品、生物活性剂、聚合物和染料的重要中间体化合物,其中伯胺是衍生化反应的关键试剂。传统合成胺的方法包括氨与烷基卤的直接反应,或通过还原硝基络合物和腈,但是这些工艺的缺陷在于产生有害的副产物、废弃物并消耗大量氢气。而利用氨气通过借氢反应机理直接胺化脂肪类醇用于合成胺被认为是一种有潜力的绿色合成路线并且已经有数千吨规模的工业化应用。然而,合成胺工业仍然面临诸多挑战,尤其是选择性合成伯胺。在此背景下,探索催化剂的本征性质与合成胺活性和选择性的依赖关系变得尤为重要。近年来,理论模拟的飞速发展使得揭示催化过程中的构效关系成为可能。本文结合第一性原理微观动力学模拟和实验,揭示了决定醇的胺化活性和选择性的关键因素并构建了其与碳和氧原子吸附能的依赖关系,实现了该反应金属合金催化剂的理性高效筛选,为实验上醇胺化催化剂的设计和升级提供了重要理论参考。

内容和讨论

A. 醇的胺化反应机理: 首先甲醇的胺化被作为模型反应进行研究,该反应机理主要包括三部分(如图1):醇的脱氢生成醛;醛与NH3迅速反应生成亚胺;亚胺的加氢生成胺。基于上述反应机理,DFT计算在九种过渡金属的密堆积表面系统展开。由于实际反应气氛中具有大量的NH3,因此催化剂的表面会覆盖一定浓度的NH3,早期的研究发现表面覆盖1/9 ML的NH3有助于准确描述胺化反应催化剂的活性顺序。其中,第二步醛与NH3反应生成亚胺被认为是气相里的快速平衡反应而且是理论模拟上的难点,该工作通过详细的反应器动力学模拟证明采用Eley–Rideal 反应机理能够准确描述该关键步骤。最终为微观动力学模拟提供了可靠的反应机理。

图1: 甲醇与氨气的胺化反应机理

B. 线性降维: 基于Norskov课题组提出的d带中心理论和scaling关系,反应过程中表面物种的吸附能并非相互独立的而是受制于scaling关系,并且基元反应的能垒也可以通过BEP关系与反应热关联,因此通过详细的描述符的筛选以及线性拟合可以实现降维和最优描述符的确认。由于该反应涉及含C、N、O的三类中间体,scaling关系可以最终将该反应中涉及的中间体的能量与C、N和O原子的吸附能进行关联。然而,包含三种不同描述符导致无法利用典型的二维火山曲线进行活性描述。深入研究发现:N原子的吸附能同样依赖于C和O原子的吸附能并且具有完美的线性关系(如图2)。而该依赖关系的物理化学根源一方面为N(3.0)原子的电负性恰好介于C(2.5)和O(3.5)原子之间,另一方面N原子的最高占据轨道能量(HOMO)可以线性分解为C和O原子的相应能量。基于上述分析,甲醇胺化反应的能量维度被降低为二(C和O原子的吸附能)。

**

图2: N原子的吸附能与C和O原子吸附能的线性关系

C. 醇的胺化活性和选择性与简单描述符的依赖关系: 基于图1的反应机理和图2的线性降维,通过构建微观动力学模型,最终建立了醇胺化的反应速率(TOF)与C和O原子吸附能的依赖关系。如图3a所示,由于对C和O原子吸附能力的区别,九种金属分布于火山图的不同位置而呈现出不同的反应活性。然而,针对醇的胺化反应,人们更加关心且关键的问题是选择性。因此,该工作进一步系统考察了导致伯胺选择性降低的不同可能性,最终发现CH3NH和CH2NH物种的耦合是导致选择性丢失的关键因素之一,并巧妙的将其与简单描述符进行了关联,最终构建了醇胺化生成伯胺的选择性与C和O原子吸附能的依赖关系如图3b。随着CH3NH/CH2NH耦合能垒(Ea)的增加,副反应变得困难,进而催化剂的选择性增加。至此,通过计算催化剂表面C和O原子的吸附能,结合图3中的火山图便可定性预测相应催化剂的醇胺化反应的活性和选择性,从而大大提高催化剂的筛选效率。

图3: 醇的胺化活性(a)和选择性(b)与C和O原子吸附能的依赖关系

D. 实验室论证预测: 基于上述活性和选择性图,通过对元素周期表中不同过渡金属进行两两排列组合,最终计算了三百余种双金属合金的C和O原子吸附能,并预测出系列有潜力的醇选择性胺化催化剂。基于上述理论预测,索尔维公司研发中心对不同合金进行了详细的实验室测试以论证图3的可靠性。具体的:本文合成了系列Co基合金催化剂,一方面由于Co金属的相对廉价,另一方面Co具有较高的活性和伯胺选择性。根据图3a的理论预测,提高Co的活性需要设法降低其O原子的吸附能,由于Pt和Pd具有弱的O吸附能力,因此将Co与Pt和Pd形成合金将可以实现反应活性提高的,然而基于图3b的信息,这将同时伴随着伯胺选择性的降低。事实上,图4和图5的实验结果也证实了CoPt和CoPd合金相比Co金属的高活性和低伯胺选择性。此外,图3表明:将Co与Ag和Ru组合将会维持甚至提高活性且提高选择性。该预测结果同样与图4和图5中的实验结果一致,并且成功设计出活性选择性都大大提高的CoAg合金催化剂。本文预测的其他有潜力的合金催化剂也将进一步推动实验上对醇胺化反应研究。

img

图4: 实验上不同催化剂醇的胺化生成伯胺的TOF

图5: 实验上不同催化剂醇的胺化生成伯胺的选择性

总结和展望

通过精确的第一性原理计算和微观动力学方法,准确的模拟了醇的胺化这一复杂反应的详细机理,揭示了决定金属合金催化剂醇的胺化活性和选择性的关键因素,并与简单易得的C和O原子吸附能进行关联,系统的实验测试论证了理论预测的可靠性。最终大大提高了醇胺化反应金属合金催化剂的筛选效率,为工业上醇胺化催化剂的设计和升级提供了参考。此外,该工作采用的改进的动力学模拟方法能够引入反应器环境对复杂反应机理的影响,一方面实现对实验条件更加精准的描述,另一方面能够准确定性预测金属合金催化剂的活性和选择性的变化趋势,为后续高通量催化剂筛选提供了基础。

作者简介

文章第一作者王涛博士目前为美国斯坦福大学SUNCAT中心博士后研究员,通讯作者Philippe Sautet为美国加州大学洛杉矶分校教授、ACS Catalysis副主编、法国科学院院士、美国能源局IMASC能源前沿研究中心副主任,长期致力于理论催化方向的基础研究。此外,法国里昂高等师范学院、法国里尔大学、比利时根特大学以及跨国公司索尔维研发中心E2P2实验室也为该工作提供了支持。本文经大师兄涛本人授权发布。

文章链接:https://www.nature.com/articles/s41929-019-0327-2

关注公众号:大师兄科研网,回复:胺化,可以获取下载链接。

工作简介

上海交通大学密西根学院朱虹老师课题组诚聘博士后 1~2 名,负责或协助团队进行金属合金体系的腐蚀性能的理论预测和腐蚀机理的研究,将与腐蚀实验团队以及世界各地的材料基因组计算研究团队开展紧密合作 (UC Berkeley, University of Michigan, UC San Diego, Georgia Tech, University of Maryland, etc),并在现有的高通量计算平台上开发针对腐蚀相关基础参数的高通量计算工作流。诚挚邀请具有理论计算、电化学背景同学的加入。

导师简介

朱虹 (hong.zhu@sjtu.edu.cn) 助理教授,上海交通大学材料基因组联合研究中心成员; 材料科学与工程学院,上海交通大学密西根联合学院双聘。

Google Scholar: https://scholar.google.com/citations?user=x1BGIfEAAAAJ

课题组主页: http://umji.sjtu.edu.cn/~hzhu

工作待遇

  1. 年收入 22w+, 享受相应上海交通大学福利待遇;
  2. 根据上海市博士后管理政策办理有关落户事宜;
  3. 优秀者可转入专职科研序列。

招聘要求

  1. 拥有材料、化学或相关专业的博士学位
  2. 具有理论模拟计算经验或具有丰富的机器学习特别是材料研究领域的经验
  3. 能够独立开展研究,在知名期刊上以第一作者发表论文
  4. 在同等条件下,优先考虑具有较高编程能力的同学

申请方式

发邮件至 hong.zhu@sjtu.edu.cn,并附上简历以及可以到站的时间。

Python现在越来越流行,连微信朋友圈都开始推广一些python培训课程广告了:俩美女一起上班,一个瞬间完成工作,另一个一脸懵逼地问道你咋这么快就干完活啦?另一个回答道:因为我学习了Python啊。是的,Python确实有这些神奇的功效,但在网络这么发达的时代,想学习python的途径有很多,参加什么割韭菜的培训班啦,找个计算机学院的男女朋友啦,向师兄师姐学习啦,自己学习琢磨啦等等。

学习Python的一个好处就是,很多懂python的都会github上分享一些自己的脚本,但这些脚本很多时候不是为你的课题量身定制的。稍微懂点可以加以修改就变成自己科研的一个利器了。更重要的好处就是广告里面所凸显的,由于计算中很多操作都是可以重复的,写个小脚本可以快速完成自己的工作任务,提高效率,省出时间该干嘛干嘛去。此外,很多复杂的任务,手动计算量很大,也得需要写点脚本来帮忙。那种瞬间完成任务的感觉,不会点程序语言的人是不会体验到的。而且,现在小学生,初中生都开始学习python了,我们也没有理由不跟上时代的潮流。

看到这里,应该能明白本文的目的。不要你交99元,也不给你优惠99元。大师兄推荐一本武林秘籍,专门针对科研人员的一本科学计算相关的Python书。里面大量例子简单易懂,解释详细,带你缓缓走入Python的程序世界,平日闲的没事加以练习,犹如武林高手给你灌输真气一般,达到计算的另一个境界,轻而易举超越那些培训班被割过的韭菜们。

获取方式,扫描右侧或者下方二维码关注公众号:大师兄科研网,回复:武林秘籍 就可以获取下载链接了。为了避免别人学会了python跟你抢男女朋友,课题组争宠,大家下载后私下练习就OK啦,有对象的,一起练习讨论效果更佳。

硕博(后)招聘:美国马萨诸塞大学洛厄尔分校车芳琳课题组

(University of Massachusetts Lowell)

课题组简介

车芳琳老师课题组刚刚成立,主要研究催化以及材料科学方面的多尺度理论模拟,在多尺度范围内(从原子到反应器尺度)研究燃料电池,电催化及新型能源、光催化相关的催化剂性质,反应机理,动力学等。研究的相关课题有:轻质烷烃的活化,CO2电催化,多相催化的线性关系以及微观动力学模拟,以及多相催化反应器中的流体动力学模拟。课题组主页:http://sites.uml.edu/fanglin_che/ (或点击左下角阅读原文)

导师简介

车老师本科毕业于大连理工大学(2008-2012),博士毕业于华盛顿州立大学(2012-2016),师从J.-S. McEwen教授。分别在加拿大多伦多E. H. Sargent课题组( 2017-2018)和美国特拉华大学Dionisios G. Vlachos课题组(2018-2019)从事博士后研究。至今以一作(含共一)在Nat. Chem. (1), Nat. communs (2), Angew. Chem. Int., Ed.(1), ACS Catal.(3), J. Catal.(1), ACS Energy Lett.(1), Adv. Mater(1), Appl. Catal. B (2)等高水平期刊发表文章若干篇。

Googl Scholar:

https://scholar.google.com/citations?hl=zh-CN&user=fiRHcIQAAAAJ

申请方式

目前长期招聘硕士,博士,以及博士后。有兴趣的可以将个人简历发送至:Fanglin_che@uml.edu

我们计算的时候,经常会遇到扩胞的需求。比如,我们要优化(4x4)的Ag(111)的一个slab。直接切一个slab拿来优化,可能会比较耗时。另一个办法就是我们先优化一个(1x1)的Ag(111)的slab,然后在将优化完的结构扩成(4x4)的,最后再优化。这样能有效地减少工作量。

扩胞的话,有很多工具可以选择,MS, P4vasp,等等可视化的软件,用鼠标点点就OK了。也可以使用一些现成的脚本,小程序。本节就介绍一个通过ASE进行扩胞的小脚本,也是本人偶然在一个网站发现的。有兴趣的可以自己看下: https://www.nsc.liu.se/~pla/blog/2013/02/26/vaspsupercells/

废话不多说,直接上例子,例子完了是脚本的具体内容。

如果你有自己的脚本或者推荐的程序,也欢迎发送到本人邮箱:lqcata@gmail.com。后面我们会逐渐扩展本节的内容。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
qli@bigbro:~/Desktop/test_expand$ ls
CONTCAR expand.py
qli@bigbro:~/Desktop/test_expand$ cat CONTCAR
Ag
1.00000000000000
3.0472766735234269 0.0000000000000000 0.0000000000000000
1.5236383367617135 2.6390190116310261 0.0000000000000000
0.0000000000000000 0.0000000000000000 22.4642729552180782
Ag
4
Selective dynamics
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 F F F
0.3333333333333357 0.3333333333333357 0.1107576902236147 F F F
0.6666666666666643 0.6666666666666643 0.2214586458624846 T T T
-0.0000000000000000 -0.0000000000000000 0.3323996157642722 T T T

0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
qli@bigbro:~/Desktop/test_expand$ python expand.py CONTCAR 2 2 1
qli@bigbro:~/Desktop/test_expand$ ls
CONTCAR expand.py POSCAR_ex
qli@bigbro:~/Desktop/test_expand$ cat POSCAR_ex
Supercell
1.0000000000000000
6.0945533470468538 0.0000000000000000 0.0000000000000000
3.0472766735234269 5.2780380232620523 0.0000000000000000
0.0000000000000000 0.0000000000000000 22.4642729552180782
16
Selective dynamics
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 F F F
0.1666666666666679 0.1666666666666679 0.1107576902236147 F F F
0.3333333333333321 0.3333333333333321 0.2214586458624846 T T T
0.0000000000000000 0.0000000000000000 0.3323996157642722 T T T
0.0000000000000000 0.5000000000000000 0.0000000000000000 F F F
0.1666666666666679 0.6666666666666679 0.1107576902236147 F F F
0.3333333333333321 0.8333333333333323 0.2214586458624846 T T T
0.0000000000000000 0.5000000000000000 0.3323996157642722 T T T
0.4999999999999999 0.0000000000000000 0.0000000000000000 F F F
0.6666666666666679 0.1666666666666679 0.1107576902236147 F F F
0.8333333333333320 0.3333333333333321 0.2214586458624846 T T T
0.4999999999999999 0.0000000000000000 0.3323996157642722 T T T
0.4999999999999999 0.5000000000000000 0.0000000000000000 F F F
0.6666666666666677 0.6666666666666679 0.1107576902236147 F F F
0.8333333333333320 0.8333333333333323 0.2214586458624846 T T T
0.4999999999999999 0.5000000000000000 0.3323996157642722 T T T

脚本内容如下:

1
2
3
4
5
6
7
8
9
#!/usr/bin/env python3
## https://www.nsc.liu.se/~pla/blog/2013/02/26/vaspsupercells/
import ase.io.vasp
import sys

cell_file = sys.argv[1]
x,y,z = [int(i) for i in sys.argv[2:5]]
cell = ase.io.vasp.read_vasp(cell_file)
ase.io.vasp.write_vasp("POSCAR_ex",cell*(x,y,z), label='Supercell',direct=True,sort=True)

使用前提是你已经安装好了ASE。

再次啰嗦一下,欢迎大家分享自己扩胞的小脚本,推荐的小程序,一起完善本节的内容。请发送到本人邮箱:lqcata@gmail.com。

img

汇集本人日常学习中的一些数据库,个人感觉不错的网站。也会不断地更新完善。如果你有好的数据库,链接,推荐的书籍,也欢迎发邮件至lqcata@gmail.com留言帮助我们更新。

结构查找

The Materials Project

https://www.materialsproject.org/

使用方法:

1) 点击右上方:Login

2) 点击: Sign in with your email address

3) 输入自己的邮箱,然后点击:send login link

4) 打开邮箱,点击链接,即可查询数据库。

ioChem-BD

http://www.iochem-bd.org/

本人老板和其他几个老板一起做的数据库,涵盖VASP,高斯等主流程序,可以上传结果。数据库依托于巴塞罗那超算中心。本人博士,博后的所有工作都在这个数据库中。

ChemSpider

http://www.chemspider.com/ RSC旗下的一个数据库,非常方便查找下载结构。

Pubchem

https://pubchem.ncbi.nlm.nih.gov/

ICSD 数据库

http://www2.fiz-karlsruhe.de/icsd_home.html , 需要账号密码

CCDC: The Cambridge Crystallographic Data Centre

https://www.ccdc.cam.ac.uk/

cccbdb 数据库

https://cccbdb.nist.gov

团簇数据库:

http://www-wales.ch.cam.ac.uk/~wales/CCD/CoCCD/cobalt.html

可能有些老了。

https://ww2.mathworks.cn/matlabcentral/fileexchange/33449-cluster-generator

Wulffman

https://www.ctcms.nist.gov/wulffman/

XPS

https://xpssimplified.com/

Cn Fullerenes

http://www.nanotube.msu.edu/fullerene/fullerene-isomers.html

热力学数据

NIST 数据库: https://physics.nist.gov/cuu/Constants/index.html

CRC Handbook of physical chemistry : http://hbcponline.com/faces/contents/ContentsSearch.xhtml 查找热力学参数,晶格常数,entropy等等。

http://www2.ucdsb.on.ca/tiss/stretton/database/inorganic_thermo.htm 无机化合物的物理和热力学相关数据

http://www.knowledgedoor.com/2/elements_handbook/cohesive_energy.html 结合能

http://www.surface-tension.de/ 表面张力

http://www.wiredchemist.com/chemistry/data/thermodynamic-data

https://atct.anl.gov/

物理常数

https://physics.nist.gov/cuu/Constants/index.html

http://wild.life.nctu.edu.tw/class/common/energy-unit-conv-table.html

http://halas.rice.edu/conversions

书籍文献下载:

google和google scholar:科学上网。

http://gen.lib.rus.ec/ 下载书籍

资源丰富,但不要贪多,更不要下一堆放到自己电脑里占地方。找几本好书,认真静下心来好好学习。

http://www.expaper.cn/ 科研速递论坛, 求助文献,书籍,成功率极高。

sci-hub 系列,现在网上各种域名,插件,桌面版的都有,我就不瞎凑热闹了。

5000人QQ群: 157099073 汇聚大家的力量,解决每一篇小文献的困难。

书籍和文献下载的原则,宁缺毋滥,认真阅读每一篇下载的文献,对自己的时间负责。别图多。

VASP 相关

VASP官网

http://www.vasp.at/, VASP最权威,最专业的学习资料都在这里面了。英文刚开始看起来有点不适应,坚持下去,慢慢就好了。不建议看除官网以外的那些乱七八糟的说明(包括我自己写的教程)

VASP 教程视频

http://www.nersc.gov/users/training/events/3-day-vasp-workshop/, A 3-day VASP Workshop at NERSC: VASP 开发者(长发大胡子那哥们)做的workshop, youtube 有视频。链接:https://www.youtube.com/playlist?list=PL20S5EeApOSumWZkzsaYxAvozjvFZ3ks4

VASP官方论坛

https://cms.mpi.univie.ac.at/vasp-forum/

VASP官方论坛,需要注册才可以提问,但一般来说你的问题已经有人提问过了,直接复制到google里面,基本都会显示该论坛的链接。可以不保存。

Henkelman 课题组

http://theory.cm.utexas.edu/henkelman/, 大量VASP相关的脚本,VTST计算过渡态,bader电荷分析,CHGCAR处理脚本等,向Henkelman致以崇高的敬意。

ASE: Atomic Simulation Environment: https://wiki.fysik.dtu.dk/ase/, 适合催化计算相关

Pymatgenhttp://pymatgen.org/, 适合材料计算的

K-pathhttps://www.materialscloud.org/work/tools/seekpath

能带Kpath选择的网站,上传自己的结构,便可以得到建议的路径。

Chemml: https://hachmannlab.github.io/chemml/

p4vasphttp://www.p4vasp.at/

RDkithttps://www.rdkit.org/

VASPkithttps://sourceforge.net/projects/vaspkit/files/ 来自中国的良心软件, VASP计算前后处理,功能强大,使用简单,老少皆宜。

VESTA http://jp-minerals.org/vesta/en/ 来自日本的良心软件

Xcrysdenhttp://www.xcrysden.org/, 用于批量做图,选择K-path做能带图,等。

Openbabelhttp://openbabel.org/wiki/Main_Page 强大的格式转换工具,它说第二,没人敢说第一。

RDkit:化学信息学重要的开源包。

程序软件相关:

Bash :http://tldp.org/HOWTO/Bash-Prog-Intro-HOWTO.html

https://learnpythonthehardway.org/

https://www.liaoxuefeng.com, 廖雪峰的网站,英文不想读的话,可以到这里学习。

https://matplotlib.org/ 用python画图

https://www.scipy.org/ python 科学运算相关的库

https://stackoverflow.com/ 一般谷歌搜索遇到出现的程序问题,都可以在这里找到答案。

https://www.geeksforgeeks.org

https://stackoverflow.com/ 一般谷歌搜索遇到出现的程序问题,都可以在这里找到答案。

大牛们的Git-Hub 网站

很多懂程序, 喜欢分享的大牛们都有自己的github网址,下面会逐渐完善。

1)https://github.com/Ionizing (www.bigbrosci.com这个网站就是这牛人帮忙搭建的)

2)https://github.com/wangvei (VASPkit的作者,王伟老师)

3)https://github.com/QijingZheng (很屌的一个大牛,虽然没见过)

4) https://github.com/obaica (只要Follow了他,基本就可以找到其他大牛的github了)

5)https://github.com/tamaswells (脚本小王子,the king of sharing)

6)https://github.com/jensengroup (总有一些人,会默默地支持穷人们的科研)

7)https://github.com/LePingKYXK

8)http://home.ustc.edu.cn/~lipai/scripts/vasp.html

理论计算,推荐书籍:

1) Density Functional Theory: A Practical Introduction

http://onlinelibrary.wiley.com/book/10.1002/9780470447710

有权限的直接下载;没有权限的不要在网上瞎求,很多公式都不全,这个要注意。

2)Atkin’s Physical Chemisty (10th)

网上到处都是,最新的是第10版,如果下载不了,看第9版即可。

3)of Modern Catalysis and Kinetics

做多相催化的强烈建议!http://onlinelibrary.wiley.com/book/10.1002/3527602658

以上书籍,如不提供下载链接,基本都可以在 http://gen.lib.rus.ec/ 这个网站获得。

溶剂相关的:

https://www.organicdivision.org/wp-content/uploads/2016/12/organic_solvents.html

http://murov.info/orgsolvents.htm

求职相关

Psik Network: http://psi-k.net/

注册后,一旦有人在网络里面发送海外职位信息,你便会收到邮件提醒。更为重要的是,网站里面有很多珍贵的学习资料。

Linkedin 领英 https://www.linkedin.com/,

有些老板也会在这上面发布招聘信息,本人当年的博士位置,就是在这里面发现的。

Research gate: https://www.researchgate.net/,

可以用来下文献,关注同行最新进展。


(七)理论计算,推荐书籍:

1 DFT: A Practical Introduction

http://onlinelibrary.wiley.com/book/10.1002/9780470447710

1)有权限的直接下载

2)没有权限的不要在网上瞎求,很多公式都不全,这个要注意。

2 Atkin’s Physical Chemisty (10th)

网上到处都是,最新的是第10版,如果下载不了,看第9版即可。

3 Concepts of Modern Catalysis and Kinetics 多相催化的强烈建议!

http://onlinelibrary.wiley.com/book/10.1002/3527602658

以上书籍,如不提供下载链接,基本都可以在 http://gen.lib.rus.ec/ 这个网站获得。


其他:

1 颜色数据库:

https://www.rapidtables.com/web/color/RGB_Color.html

2 元素周期表:

http://www.rsc.org/periodic-table

https://www.ptable.com/

http://www.periodictable.com/

3 Latex相关:

本人的文章,除合作的外, 基本都是用latex写的,但是推荐网址的话,确实想不起有什么来。找一个latex入门pdf,练上几天,遇到问题google搜索,慢慢就好了。如果非要推荐一个网址的话,那么就选这个: https://tex.stackexchange.com/ 和 www.google.com


img

img

也欢迎大家关注微信公众号

前面我们介绍完了算表面吸附,以及过渡态的一些基本的操作,和注意事项。当我们面对一个新的课题时,就需要运用所学到的这些技能来完成所必需的计算,来验证我们的想法,思路等。后面几节,主要参考本人2014年发表的一篇关于甲醇分解反应计算的文章来介绍一下怎么运用所学习到的这些本领。(https://pubs.acs.org/doi/abs/10.1021/cs501698w)

搬砖这个词用来形容我们完成一个课题的计算过程,简直不能再形象不过了。

首先,搬砖这一过程, 我们要明确几个要素:

1)老板:让你搬砖的人

2)搬砖工人:你自己

3)任务:搬砖

4)目的:盖房子

下面,你要明确一点:不会计算的搬砖工不是好搬砖工。

这里说的计算,不是我们的计算化学中的计算。而是把我们的工作如何一步一步分解,结合自己的实际情况,算算需要多少的时间完活。毕竟房子总有盖完的时候,搬砖工不能一直守着。所以,当老板告诉你要盖一个大房子的时候,需要多少砖这一点你要清楚。然而,很多人都是第一次搬砖,或者之前没搬过表面计算的砖,对上面说的这一点并没有一个很清晰的概念。所以,摸清自己的现实情况,这也就是本节的重点啦。

人有高低胖瘦之分,搬砖工也不例外,有的劲大,有的劲小。劲大搬的砖多,劲小搬的少。最后到底能搬多少砖,一方面取决于你的身体条件(自身因素,搬砖前的学习经历),另一方面取决于你的工资伙食情况(老板因素)。身体倍儿棒,老板钱儿多,这样的情况通常不会看本教程,有啥问题自己(组里)就能解决。身体倍儿棒,老板没钱,或者身体倍儿不棒,老板有钱的。就需要认真动动脑子,分析下当前的形势。没劲儿又没钱的,就应该更加注意了,需要更多地在脑力上下功夫。

比如,这篇文章中,老板让我们算一遍甲醇在四个金属上的分解反应(倒过来就是合成反应)。

第一步:分析下需要搬多少砖块

别看甲醇分子简单,要彻底分解成基本的C、O、H,中间有很多的基元反应要计算。从上面的图可以发现:有三类的断裂反应:C-H,O-H,C-O键的断裂,第一步这三种都可能发生的。这一步的产物在第二步中又可以发生哪一类的反应,依次类推。 最终我们可以估算一下需要计算的基元反应,以及中间体的结构。 所以,计算之前,多分析下这些可能的过程,基本的框架,多少反应,多少物种要有个概念。可以自己尝试着画,亦可以参考文献中寻找答案。

2)我们取的四个金属:Cu,Ru,Pt,Pd稳定的表面。所以前面的计算要乘以4。

3)前面说的是最理想的情况,而实际情况则是:

A) 对每一个表面吸附的物种来说,我们要尝试不同的位点,来获取一个稳定的构型;

B)对基元反应来说,这些过渡态,100%不可能都100%一次性找到,还要考虑不收敛,不同的可能性等;

这些都会使得计算量增加。所以,第一步,大体上有多少需要计算的东西,应该有个框架。

第二步:确定一个合理的搬砖方案。

计算资源是老板提供的,也就是你的伙食。伙食好,干活就有劲。但不管有没有劲,都不愿意大热天地一个劲地搬砖。所以这一步,我们要充分结合自己现有的计算资源,来制定一个合理的策略,用最少地劲搬最多的转。该怎么做呢?

A) 善于利用已经发表的文章的数据,比如,有些结构可以在支持信息里面找到,还有一些数据库里可以下载,也可以问作者要(可能性比较低), 可以理解为找朋友一起帮忙搬砖;

B) 选择合适的slab模型:3x3还是4x4,slab取四层,还是取5层?这些是课题开始之前一些基本的测试工作,可以参考文献中别人的做法,也可以根据测试的结果自己合理选择。可以理解为:搬多大的砖块。

C) 选择合适的参数: 计算参数不对,很可能导致计算的结构或者能量有问题。 但这些都要具体分析,有些能量有问题,但结构还算OK的可以调整下参数,作为一个理想的初始结构继续用。可以理解为:半路翻车,捡起来那些没摔碎的转,继续搬。

D)善于使用前面我们介绍的快速获取理想初始结构的方法。可以理解为先用车把砖块搬到离工地最近的地方,省去往返来回跑的劲。

上面说的这些,尽可能在课题完全开展前做到位,因为它们不会花费很大的劲去做,但会节省后面很大的劲。而且,伙食好坏(计算资源给不给力),测试的过程一目了然。

第三步:结合自己的体力,伙食,认真搬砖。

这一步就简单啦,体力好同时进行,左手生擒中间体,右手活拿过渡态。体力不好,俩手搬一块砖,累了饿了(没资源)就一边凉快去。

小结

不过对于新手来说,课题刚刚开始就真正掌握上面的这几点,难度有些大,所以建议先做些前提的准备工作:

1) 熟悉自己的服务器,运算能力;

2)多查找文献,整理需要计算的框架,幸运的话,可以从支持信息直接找到结构;

3)多用小体系做测试,测试完了要对结果多思考总结。不要上来就狂交任务,最后把服务器累个半死,还不出什么好的结果。

伴随着Python2的落幕,Python3的兴起,很多基于2的软件也面临着着要么被遗忘,要么主动拥抱Python3的选择。其中就有我们最喜欢的p4vasp。https://github.com/orest-d/p4vasp 或者www.p4vasp.at。本人不认识作者,也不知道啥时候出新的版本。使用p4vasp需要有python2,但其他软件又需要python3。为此我们不得不做出一些选择, 最残忍的莫过于弃用p4vasp啦。本节介绍两个办法:

1) 修改p4vasp运行默认的python版本

2)通过update-alternatives切换系统默认的python版本。

办法-1

这个方法是公众号的牛牛(刘科)留言提出来的:

第一步: 修改p4v脚本:将最后一行中的python 改成python2

1
2
3
4
5
6
7
qli@lhz:~$ sudo gedit /usr/bin/p4v 
qli@lhz:~$ cat /usr/bin/p4v
#!/bin/sh
export LD_PRELOAD=libstdc++.so.6
export UBUNTU_MENUPROXY=0
export P4VASP_HOME=/usr/share/p4vasp
exec /usr/bin/python2 /usr/share/p4vasp/p4v.py "$@"

第二步 修改p4v.py, 将第一行中的python 改成python2

1
2
3
4
qli@lhz:~$ sudo gedit /usr/share/p4vasp/p4v.py
qli@lhz:~$ head -n 2 /usr/share/p4vasp/p4v.py
#!/usr/bin/env python2
# Copyright (C) 2003 Orest Dubay <orest.dubay@univie.ac.at>

运行实例:

1
2
3
4
5
6
7
qli@lhz:~$ python  --version 
Python 3.7.3
qli@lhz:~$ p4v
Gtk-Message: 08:12:51.061: Failed to load module "canberra-gtk-module"
(p4v.py:27951): libglade-WARNING **: 08:12:51.176: Could not load support for `gnome': libgnome.so: cannot open shared object file: No such file or directory
(p4v.py:27951): libglade-WARNING **: 08:12:51.519: unknown attribute `swapped' for <signal>.
Exception TypeError: "'NoneType' object is not callable" in <object repr() failed> ignored

办法-2

百度或者google里面搜一下:update-alternatives python 就会得到一系列的操作教程。本文不再瞎扯。

第一步

查看自己系统中python2和python3的版本,终端里面输入:python然后摁tab键,显示如下,本人安装了python2.7 和python3.7。

1
2
3
qli@lhz:~$ 
qli@lhz:~$ python
python python2.7-config python3 python3.7m python3-jsonschema python3-pbr pythontex python2 python2-config python3.7 python3.7m-config python3m python3-unit2 pythontex3 python2.7 python2-jsonschema python3.7-config python3-config python3m-config python-config

第二步

设置update-alternatives:在终端里面,依次输入下面两个命令行:

1
qli@lhz:~$ sudo update-alternatives --install /usr/bin/python python /usr/bin/python2.7 1
1
qli@lhz:~$ sudo update-alternatives --install /usr/bin/python python /usr/bin/python3.7 2

这样就设置OK啦。更加具体的,网上已经泛滥了

第三步

切换python版本,使用下面的命令:

1
qli@lhz:~$ sudo update-alternatives --config python

下面是具体运行的例子,首先本人将上面的命令行保存到.bashrc文件中。

1
2
3
4
5
6
7
8
9
10
11
12
13
qli@lhz:~$ grep pversion .bashrc
alias pversion='sudo update-alternatives --config python'
qli@lhz:~$
qli@lhz:~$ python --version
Python 2.7.16
qli@lhz:~$ p4v
Gtk-Message: 08:03:38.110: Failed to load module "canberra-gtk-module"
(p4v.py:3156): libglade-WARNING **: 08:03:38.346: Could not load support for `gnome': libgnome.so: cannot open shared object file: No such file or directory
(p4v.py:3156): libglade-WARNING **: 08:03:38.782: unknown attribute `swapped' for <signal>.
Exception TypeError: "'NoneType' object is not callable" in <object repr() failed> ignored
Exception TypeError: "'NoneType' object is not callable" in <object repr() failed> ignored
Exception TypeError: "'NoneType' object is not callable" in <object repr() failed> ignored
Exception TypeError: "'NoneType' object is not callable" in <object repr() failed> ignored

当前版本是python2.7,可以正常运行p4vasp。

1
2
3
4
5
6
7
8
9
10
11
12
13
qli@lhz:~$ pversion 
[sudo] password for qli:
There are 2 choices for the alternative python (providing /usr/bin/python).

Selection Path Priority Status
------------------------------------------------------------
0 /usr/bin/python2.7 2 auto mode
* 1 /usr/bin/python2.7 2 manual mode
2 /usr/bin/python3.7 1 manual mode

Press <enter> to keep the current choice[*], or type selection number: 2
update-alternatives: using /usr/bin/python3.7 to provide /usr/bin/python (python) in manual mode

当前的python版本,会在最左侧那一列用* 标出来。换成python3的话,输入 左侧所对应的数值就OK了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
qli@lhz:~$ python --version 
Python 3.7.3
qli@lhz:~$ p4v
File "/usr/share/p4vasp/p4v.py", line 9
can get source code from http://www.pygtk.org """
^
SyntaxError: invalid syntax
qli@lhz:~$ pversion
There are 2 choices for the alternative python (providing /usr/bin/python).

Selection Path Priority Status
------------------------------------------------------------
0 /usr/bin/python2.7 2 auto mode
1 /usr/bin/python2.7 2 manual mode
* 2 /usr/bin/python3.7 1 manual mode

Press <enter> to keep the current choice[*], or type selection number: 1
update-alternatives: using /usr/bin/python2.7 to provide /usr/bin/python (python) in manual mode

顺利将python切换到版本3,但p4vasp这时候就不能用了。

小结

希望p4vasp赶紧出基于python3的版本,如果不出的话,大家就先这样凑活着用吧。如果牛牛们还有更好的解决办法,可以将教程稍微整理下,发送到邮箱lqcata@gmail.com,本人会继续更新到本节当中。

前面一节,我们讲解了对于扩散这一类反应中,nebmake.pl这个命令生成IMAGES时候的一个坑。通过学习,要了解到一些脚本或者程序本身存在一些缺陷,我们在使用的时候,要避免盲目相信,直接把脚本的结构直接拿来用。总之,原则就是尽可能获取一些具有明确物理化学意义的,比较理想的初始构型。本节,讲解一下旋转过程中,nebmake.pl的坑。这里与其说是坑,不如说是nebmake.pl不适用的情况。因为涉及到旋转过程的时候,一般得到的IMAGES的结构都不咋地,需要自己认真检查,微调下结构或者重新搭结构。

这里我们讲一个极端的例子:乙烷分子的旋转

从一个交叉式的构象到另一个交叉式的构象,两个甲基绕C—C键旋转120°。如果每隔15°插一个点, 正常的结果应该如下图所示:

但是,当使用nebmake.pl这个方法产生的IMAGES结构如下图:

可以看到,

1)初始结构中,并没有旋转的效果,H原子走的是直线的路径。

2)IMAGE中,H和C原子的距离非常小,仅仅贴在一起了,显然这样的结构非常不理想。

如果使用这样的结构进行计算,第一个离子步结束后,计算出来原子间的作用力很强,会导致后面计算中分子直接散架,而这些散架的结构通常都不收敛,如果不及时检查结构,及时杀死,它会在服务器上一直就这样算着,而你还在傻傻地啃着西瓜,聊着QQ,等待结果。

解决办法:

知道有这个坑之后,怎么躲就好办多了。

1)对于类似的旋转结构,自己手动搭建;

2)使用其他的高级点的生成IMAGES的方法,例如IDPP。(https://wiki.fysik.dtu.dk/ase/tutorials/neb/idpp.html)

3)自己写脚本实现旋转的过程。

小节:

这两节简单介绍了一下生成NEB计算IMAGES需要注意的地方,不管咋地,原则还是要再啰嗦一遍:尽可能得到具有理想物理化学意义的初始结构。毕竟好的开始是成功的一半。有兴趣的可以算一下这两节例子,加深一下自己的印象。有大佬公众号留言说贴自己的代码教程怒怼IDPP,希望看到的可以联系俺(lqcata@gmail.com)。

<<<<<<< HEAD
最近整理材料,发现遗落在角落里的一个python脚本(cssm.py),可以批量切常见金属稳定表面的slab模型。 cssm.py 是:cleave_stable_surfaces_from_metals的缩写。废话不多说,下面的是在天河II上运行的实例。

阅读全文 »

被批评公众号长草了,简单介绍一下ASE转换文件的小技巧,顺便锄下草。顺利学会下面操作的前提:

1) 已经在自己电脑上安装好了ASE。

2) 电脑里面有POSCARtoolkit.py脚本。见(https://www.bigbrosci.com/2018/11/11/ex73/)

终端操作

下面是终端里面转换格式的懒人命令:

1
ase gui  file_from  -o POSCAR 

1 ase gui 是我们转化文件格式的命令。

2 file_from 是你想要转换的文件,可以为cif, xyz, xsd等。ASE支持很多格式文件的读取和输出,具体的大家可以参考下面网址:

a) https://wiki.fysik.dtu.dk/ase/ase/io/io.html#module-ase.io

b)https://wiki.fysik.dtu.dk/ase/_modules/ase/io/formats.html

3 -o 代表输出文件,这里我们写的是POSCAR,也可以写成CONTCAR。

如果:

1)想转换成xyz文件,改成 -o XXX.xyz;

2)想转换成cif文件,改成 -o XXX.cif;

总之,

想转什么文件就写什么文件;

想转啥格式的就写啥后缀。

如果失败的话,可能不支持,那就接着用自己的老方法吧。

个人使用经验:

本人的计算模型经常在xyzPOSCAR两个格式之间切换。下面分析下ASE转换的优点和缺点。

1)从POSCAR到xyz文件(优点):

在xyz文件的第二行中,ASE会写入一些体系的周期性的信息,这一点非常棒。方便从xyz再转换回POSCAR或者CONTCAR。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
qli@bigbro:~/Desktop/test$ ls
CONTCAR
qli@bigbro:~/Desktop/test$ cat CONTCAR
Au\(1\1\1)
1.00000000000000
20.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 21.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 22.0000000000000000
C H
2 4
Selective dynamics
Direct
0.1615221012872249 0.4048804874587976 0.3815091574487524 T T T
0.1387568880386591 0.3452706062587786 0.3846145011058063 T T T
0.2143704422309978 0.4168857617494277 0.3865629189119751 T T T
0.1712711104214057 0.3044252978921785 0.3929860718019034 T T T
0.1271257560199507 0.4443139925926637 0.3733323288922066 T T T
0.0855583096767520 0.3352469534052885 0.3790158532393584 T T T

qli@bigbro:~/Desktop/test$ ls
CONTCAR
qli@bigbro:~/Desktop/test$ ase gui CONTCAR -o c2h4.xyz
/home/qli/anaconda2/lib/python2.7/site-packages/ase/gui/ag.py:86: UserWarning: You should be using "ase convert ..." instead!
warnings.warn('You should be using "ase convert ..." instead!')
qli@bigbro:~/Desktop/test$ cat c2h4.xyz
6
Lattice="20.0 0.0 0.0 0.0 21.0 0.0 0.0 0.0 22.0" Properties=species:S:1:pos:R:3 pbc="T T T"
C 3.23044203 8.50249024 8.39320146
C 2.77513776 7.25068273 8.46151902
H 4.28740884 8.75460100 8.50438422
H 3.42542221 6.39293126 8.64569358
H 2.54251512 9.33059384 8.21331124
H 1.71116619 7.04018602 8.33834877

Lattice=XXX 这一行可以帮助你从xyz转回POSCAR。另外,转换的时候,会出现一个warning信息,装看不见就行了。

2)从xyz到POSCAR

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
qli@bigbro:~/Desktop/test$ ase gui c2h4.xyz -o POSCAR
/home/qli/anaconda2/lib/python2.7/site-packages/ase/gui/ag.py:86: UserWarning: You should be using "ase convert ..." instead!
warnings.warn('You should be using "ase convert ..." instead!')
qli@bigbro:~/Desktop/test$ cat POSCAR -n
1 C H
2 1.0000000000000000
3 20.0000000000000000 0.0000000000000000 0.0000000000000000
4 0.0000000000000000 21.0000000000000000 0.0000000000000000
5 0.0000000000000000 0.0000000000000000 22.0000000000000000
6 2 4
7 Cartesian
8 3.2304420299999999 8.5024902400000002 8.3932014600000002
9 2.7751377599999998 7.2506827300000003 8.4615190200000008
10 4.2874088400000003 8.7546009999999992 8.5043842200000004
11 3.4254222099999998 6.3929312600000001 8.6456935799999997
12 2.5425151200000000 9.3305938400000006 8.2133112399999995
13 1.7111661899999999 7.0401860200000002 8.3383487699999996

从xyz转化成POSCAR的时候,ASE的优点:

1)按照Cartesian格式输出。

2)元素行在第一行。

缺点:

1)输出文件还是老式的VASP格式(vasp 4.6)。第6行前面少了元素哪一行。

解决办法:ASE把这一行写到第一行了。我们直接在vim里面,把第一行复制到第6行前面就OK。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 1  C  H
2 1.0000000000000000
3 20.0000000000000000 0.0000000000000000 0.0000000000000000
4 0.0000000000000000 21.0000000000000000 0.0000000000000000
5 0.0000000000000000 0.0000000000000000 22.0000000000000000
6 C H
7 2 4
8 Cartesian
9 3.2304420299999999 8.5024902400000002 8.3932014600000002
10 2.7751377599999998 7.2506827300000003 8.4615190200000008
11 4.2874088400000003 8.7546009999999992 8.5043842200000004
12 3.4254222099999998 6.3929312600000001 8.6456935799999997
13 2.5425151200000000 9.3305938400000006 8.2133112399999995
14 1.7111661899999999 7.0401860200000002 8.3383487699999996

缺点2)少了固定元素的哪一行,以及固定坐标的’T’ 和 ‘F’ 。

解决办法:通过前面介绍的POSCARtoolkit.py脚本来解决,运行下面的命令,输入固定的层数即可。具体用法见(https://www.bigbrosci.com/2018/11/11/ex73/)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
qli@bigbro:~/Desktop/test$ ls
c2h4.xyz CONTCAR POSCAR
qli@bigbro:~/Desktop/test$ POSCARtoolkit.py -i POSCAR -f
x
x
x
-----------------------------------------------------------
Cartesian Coordinates found, only for fixing atoms!
Then type how many layers to be fixed, from bottom to top.
-----------------------------------------------------------

Found 1 layers, choose how many layers to be fixed------> 0

-----------------------------------------------------------
POSCAR with Cartesian Coordiations is named as POSCAR_C
x
x
x
qli@bigbro:~/Desktop/test$ ls
c2h4.xyz CONTCAR POSCAR POSCAR_C
qli@bigbro:~/Desktop/test$ cat POSCAR_C
C H
1.0000000000000000
20.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 21.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 22.0000000000000000
C H
2 4
Selective
Cartesian
+3.2304420300 +8.5024902400 +8.3932014600 T T T
+2.7751377600 +7.2506827300 +8.4615190200 T T T
+4.2874088400 +8.7546010000 +8.5043842200 T T T
x
x
x

小结:

个人还是比较喜欢用ASE来进行文件格式的转换,方便实用。方法千万种,能找到适合自己的就行。此外,尝试新的方法的时候,一定要认真研究分析下输出的文件(脚本、程序或多或少都会有些bug),保证对于自己的体系不出错后再继续使用,安全第一。此外,如果你有自己独特的、简单粗暴的、傻瓜式的、不通过鼠标敲敲点点的转换方法,也欢迎分享。

算过渡态,光知道检查虚频是不对的。如果算出来虚频多,很大程度是因为你的NEB初始结构有问题,也就是说你的过渡态路径背后的物理化学意义不是那么地理想。而对于NEB的初始结构,大部分人都是通过VTST的nebmake.pl脚本来实现的。

1
nebmake.pl IS FS 8 

分析上面这个命令:有4个关键的信息:

  • 脚本:nebmake.pl
  • 初始结构: IS
  • 末态结构: FS
  • 插的点数:8

本节默认你优化好了IS,FS,并且插8个点,分析一下脚本nebmake.pl的2个坑。其实,在Density Functional Theory: A Practical Introduction 这本书的第6章,其中的一个坑已经讲过了。强烈建议新手,老手认真看看这一章。

想知道这两个坑,首先我们先分析下nebmake.pl的工作原理。简举个化的一维例子:一个线段两端的坐标是x1和x2,我们把这个线段分成n份,并获得每一段的起始点的坐标。想必大家都知道怎么弄。

x_i = x0 + i *(x2-x1)/(n+1)

如果扩展到三维的xyz坐标,分别对y、z进行同样的处理。我们就得到了初始结构和末态结构之间这些的IMAGES的坐标。具体见:第六章的148页。

知道了原理,现在就可以分析其中的两个坑了。

坑1

前面讲的原子在表面扩散的例子,如果原子从fcc扩散到hcp。fcc和hcp的结构中,原子在z方向的坐标基本是相同的。如果运行nebmake.pl 这个命令,就会出现这样的一种情况,所有的IMAGE结构中,z的坐标几乎不变,也就是原子在表面横着走。而实际情况呢?原子在扩散的过程中,z方向的坐标也发生相应的变化。就如同你从山的一边爬到另一边,虽然海拔基本没变,但爬山这个过程,有升,下山过程,有降。所以,这种情况下,直接用nebmake.pl插的点在z方向上的物理意义并不准确。稍微扩展下,在x或者y方向上也可能会发生类似的情况,如果IS和FS在某一维度的变化很小时,一定要注意这一维度上的物理意义是否可以被表现出来。

既然知道了这个坑,我们该怎么填呢?

方法1:

管它是不是坑,或者不知道这个坑,直接开车冲过去。这是很多新手常见的做法,虽然有时候可以开过去,但不推荐,毕竟也会溅一车泥,搞不好坑大了还会掉进去。

方法2:

学习书中的例子,在计算之前先手动修改下结构。比如插了8个点,把第2和第7个IMAGES中的坐标向上移动0.1$\AA$。第3和6个向上移动0.15$\AA$,第4和5移动0.2$\AA$。这样做的好处就是,初始结构在z方向上具有更好的物理意义,使得计算收敛更快。

上图是书中的一个例子,最上面的曲线是没有调节z坐标时,初始结构所对应的能量;中间的为调节初始结构的z坐标后,初始结构的能量;最下面的是NEB计算完成之后,各个IMAGE的能量。从能量上也可以看出来,如果z坐标我们不修改的话,体系能量与稳定的相差甚远,间接告诉我们可能需要更多的优化步数来收敛。这个强烈建议大家自己亲手算一算。体系简单,不耗费那么多机时,有助于加深对NEB的认识和学习。

前面几节,我们从gamma点到3x3x1计算过渡态。首先复习一下计算的流程:

1) 我们计算的是金属表面上H原子的扩散;

2) 使用gamma点计算的时候,slab固定住了;

3)gamma点计算结束后,检查结果:

  • 看到有一个漂亮的NEB图;
  • 能量变化也是稳稳妥妥滴;
  • 结构查看也没啥异常情况;

4) 在确认第3)步检查OK之后,将gamma点的计算备份;增大K点至3x3x1继续算。

5)计算结束后,重复第3)步的检查,确认没啥问题。

继续算

前面这么做很啰嗦(大师兄本人也很啰嗦),目的只有一个:用最少的机时获取最好的NEB初始结构。当我们完成这一步之后,就可以再继续下面的操作:

  • 备份3x3x1的计算;

  • 增大K点至5x5x1;

  • 放开表面的原子;(POSCARtoolkit.py
  • 继续算,结果如下:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ ls
00 02 04 06 08 INCAR POTCAR job_sub movie neb.dat slurm-1133307.out spline.dat vasprun.xml 01 03 05 07 09 KPOINTS exts.dat mep.eps movie.vasp nebef.dat slurm-1133315.out vaspgr
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ cat KPOINTS
K-POINTS
0
Gamma
5 5 1
0 0 0
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ tail 02/OUTCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ tail 02/OUTCAR
User time (sec): 3864.268
System time (sec): 15.316
Elapsed time (sec): 3905.371

Maximum memory used (kb): 466776.
Average memory used (kb): 0.

Minor page faults: 816554
Major page faults: 29
Voluntary context switches: 29864
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ tail 02/OSZICAR -n 1
10 F= -.32111627E+03 E0= -.32109808E+03 d E =-.130617E-03
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ ta.sh
00 -321.15187450
01 -321.13494327
02 -321.09807672
03 -321.06314408
04 -321.03855868
05 -321.02460382
06 -321.01884134
07 -321.02967204
08 -321.06853536
09 -321.09850728

从上面可以看到,顶点在06的位置,差不多就是过渡态了。

验证过渡态

那么我们算出来的过渡态到底对不对呢?

上面,我们通过能量分析,06这个Image就是过渡态了;但06就真的是过渡态吗?下面我们需要做2件事情:

1) 查看结构:

为了区分明显,上图中桥式位置两端的Ru原子,用暗红色标记出来。可以看出,06结构中,H原子在桥式的位置上;是我们想要的过渡态。结构这一关也过了。

2)频率分析:对于一个基元反应的过渡态来说,会有一个对应的虚频。因此,我们将06的CONTCAR单独取出来,做一个频率分析:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79$ mkdir freq && cp 06/CONTCAR freq/POSCAR && cd freq &&sed -i '10,27s/T/F/g' POSCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ cp ../INCAR .
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ vi INCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ kpoints.sh 1 1 1
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ pospot.sh
Generating NEW POTCAR...
************************
Done
************************
NEW POTCAR containes....
Ru H
************************
Elements in POSCAR
************************
Ru H
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ qsub
Submitted batch job 1179250
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ ls
CHG CONTCAR DYNMAT IBZKPT KPOINTS OUTCAR POSCAR REPORT XDATCAR slurm-1179250.out CHGCAR DOSCAR EIGENVAL INCAR OSZICAR PCDAT POTCAR WAVECAR job_sub vasprun.xml
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$ grep cm-1 OUTCAR
1 f = 36.804765 THz 231.251161 2PiTHz 1227.674788 cm-1 152.212329 meV
2 f = 30.534675 THz 191.855022 2PiTHz 1018.527097 cm-1 126.281311 meV
3 f/i= 9.822527 THz 61.716759 2PiTHz 327.644231 cm-1 40.622722 meV
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex79/freq$
  • 创建freq文件夹,将06的CONTCAR复制到freq中;
  • 将表面的Ru原子固定住;
  • 复制NEB计算的INCAR过来,修改频率计算相应的参数;
  • 使用gamma点算频率;
  • 生成POTCAR;
  • 提交任务。
  • 查看结果,有一个虚频。计算完事!

注意:

  • 不会算的看前面频率计算的内容;

  • 自己测试下不同K点:2x2x1;3x3x1; 5x5x1计算出来的零点能有什么区别?频率大小有什么区别。

小结:

拖拖拉拉,磨磨唧唧,终于把H原子在表面过渡态计算的NEB讲的差不多了。如果你一路沿着教程走过来,最终看到一个虚频对的结果时,会发自内心由衷地一笑:啊!啊!啊!过渡态原来这么简单。为避免这种情况的发生,我先给你脑勺拍一板凳,后面事还多着呢。还是那句话:做任何一个事情,从你认为它简单的那一刻起,你就输了。算过渡态,

1)要把初始结构,末态结构优化好;

2)合理采用粗糙的模型,来检查自己设想的反应路径;

  • 粗糙的模型:一方面指的是体系的大小,一方面是计算参数的设置;
  • 切忌直接上来硬算,要不然把服务器累个半死,却得不到多少好的结果。

3) 从能量,结构,虚频等多个角度去分析你的结构。

本节计算的文件已经打包上传到百度网盘,由于版权,压缩包里面的POTCAR就不放了,大家自己生成一下。下载链接:链接:https://pan.baidu.com/s/1sfOLxTB5Rdr5Il7vNCOkyA 提取码:m591

什么是ASE?

ASE 是Atomic Simulation Environment的简称。详细的介绍大家可以浏览官网:https://wiki.fysik.dtu.dk/ase/about.html。ASE是一个基于python的程序库,可以做的事情很多,不同程序任务的设置,结构的搭建,结构分析都可以做。

为什么用ASE?

  • ASE里面包含了很多的功能,其中个人认为最为方便地就是计算模型不同文件类型的转换。比如将cif文件转化为POSCAR,将POSCAR转化为xyz文件,等等;

  • 可以在学习ASE的时候,顺便学习python的一些知识;

  • 使用ASE分析计算结果;

  • 使用ASE学习做计算;这里分享江山推荐的一本书:https://github.com/jkitchin/dft-book

怎么安装ASE?

如果你用的是国防科大吕梁超算中心,管理员已经安装好了ASE,那么需要你做的只有一步,添加ASE的路径到.bashrc文件中就可以直接用了。

1
2
3
4
5
6
iciq-lq@ln3:/THFS/opt/python3.6/bin$ ls ase*
ase ase-build ase-db ase-gui ase-info ase-run
iciq-lq@ln3:/THFS/opt/python3.6/bin$ pwd
/THFS/opt/python3.6/bin
iciq-lq@ln3:/THFS/opt/python3.6/bin$ tail -n 1 ~/.bashrc
export PATH=$PATH:/THFS/opt/python3.6/bin

如果你想在自己电脑上安装ASE。安装过程很简单,详细的信息请参考:https://wiki.fysik.dtu.dk/ase/install.html。Mac太贵,没尝试过,Windows系统本人不太熟悉,下面附上在Windows10中Ubuntu寄生系统的安装过程(跟正儿八经Ubuntu系统的安装流程一样); 总共就三步。

  • 安装pip:如果已经安装了pip,则就剩2步了。

    1
    2
    sudo apt install python-pip
    sudo apt install python3-pip

    一个是基于python2,一个是python3的。

  • 安装ASE:

    1
    pip install --upgrade --user ase

  • 添加ASE的路径到.bashrc文件。

1
export PATH=$PATH:/home/bigbro/.local/bin

在Terminal里面输入ase,然后摁下tab键:如果跟下面的输出一样,说明八九不离十了。

1
2
3
4
bigbro@DESKTOP-S6R3TUP:~$ ase
ase ase-build ase-db ase-gui ase-info ase-run ase_figures.sh aselite.py aselite.pyc
bigbro@DESKTOP-S6R3TUP:~$ ase --version
ase-3.17.0

前面一节我们通过粗糙的方法演示了一遍怎么在超算上跑过渡态的NEB计算。一般来说,通过这一步我们可以大体上确定这个基元反应中的路径是不是靠谱的了。不论你得到一个粗糙的NEB能量曲线:

亦或者一个比较精细的:

这都表明,你快要得到这个过渡态结构了。在进行后面操作讲解之前,我们先做2个事情:

1) 由于本人也是老司机了,算过渡态的时候遇到的那些零零碎碎的问题一时也很难记起来,如果你正在算过渡态,并且遇到了问题和错误,可以将计算打包,以及遇到的问题发给我。如果本人有幸可以帮你解决问题,那么,作为交换条件,这个例子就写在后续的教程里面。有下面几点需要注意的:

  • i)请使用邮件发送例子(lqcata@gmail.com),尽可能保证计算的文件完整,千万不要删掉某些文件后再打包发给我(CHG,CHGCAR,WAVECAR除外)。如果例子太大,可以上传到百度网盘,然后将链接和你遇到的问题发给我。

  • ii)本人最近很忙,可能回复有些慢,请见谅。如果你看到我在QQ群里面瞎BB的话,可以提醒下我去查看邮件。

  • iii)只想解决自己问题,不愿意分享自己例子的,请不要发邮件给我。

2)第二件事就是分享一个自己平时常用的脚本,这个脚本主要是将前面的计算步骤保存起来,然后将CONTCAR复制成POSCAR用以下一步的计算。脚本的名字为:save_calculations.sh.

1
2
3
4
5
6
7
8
9
10
#!/usr/bin/env bash 

mv POSCAR POSCAR-$1
mv OUTCAR OUTCAR-$1
mv OSZICAR OSZICAR-$1
mv vasprun.xml vasprun.xml-$1
mv EIGENVAL EIGENVAL-$1
mv IBZKPT IBZKPT-$1
cp CONTCAR CONTCAR-$1
mv CONTCAR POSCAR

脚本演示:

1
2
3
4
5
6
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/test$ ls
CONTCAR DOSCAR DYNMAT IBZKPT INCAR KPOINTS OSZICAR OUTCAR p4vasp.log POSCAR POTCAR sub12 sub24 vasprun.xml
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/test$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/test$ save_calculations.sh gamma
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/test$ ls
DOSCAR DYNMAT IBZKPT INCAR KPOINTS OSZICAR-gamma OUTCAR-gamma p4vasp.log POSCAR POSCAR-gamma POTCAR sub12 sub24 vasprun.xml-gamma

脚本说明:

  • 脚本内容很粗暴直接,大家嫌不好看也可以结合下for循环改进地更加简洁些,本人也懒得改,能实现目的就OK了;
  • 效果很简单,就是把之前计算的一些结果重新命名一下,大家可以根据自己的爱好或者习惯,将gamma换成其他的。

提交任务

我们把参数设置地稍微精细一下,然后再跑一遍NEB。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ ls
00 01 02 03 04 05 06 07 08 09 INCAR job_sub KPOINTS NEB.pdb POTCAR slurm-995085.out vasprun.xml
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ ls *
INCAR job_sub KPOINTS NEB.pdb POTCAR slurm-995085.out vasprun.xml

00:
POSCAR

01:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT WAVECAR XDATCAR
.
.
.
08:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

09:
POSCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ for i in 0{1..8}; do cd $i ; save_calculations.sh gamma ; cd - ; done
mv: cannot stat `vasprun.xml': No such file or directory
/THFS/home/iciq-lq/LVASPTHW/ex78
.
.
.
mv: cannot stat `vasprun.xml': No such file or directory
/THFS/home/iciq-lq/LVASPTHW/ex78
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ ls *
INCAR job_sub KPOINTS NEB.pdb POTCAR slurm-995085.out vasprun.xml

00:
POSCAR

01:
CHG CHGCAR CONTCAR-gamma DOSCAR EIGENVAL-gamma IBZKPT-gamma OSZICAR-gamma OUTCAR-gamma PCDAT POSCAR POSCAR-gamma REPORT WAVECAR XDATCAR
.
.
.
08:
CHG CHGCAR CONTCAR-gamma DOSCAR EIGENVAL-gamma IBZKPT-gamma OSZICAR-gamma OUTCAR-gamma PCDAT POSCAR POSCAR-gamma REPORT stdout WAVECAR XDATCAR

09:
POSCAR

iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ kpoints.sh 3 3 1
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ yhbatch -p gsc -N 2 -J test job_sub
Submitted batch job 1133307
  • 这里我们把01到08中的一些输出文件都保存了一下,然后将CONTCAR复制成POSCAR用来进行下一步的计算;
  • 运行的时候,会报错mv: cannot stat ‘vasprun.xml’: No such file or directory,不用管就行,因为images的文件夹中没有vasprun.xml`文件,嫌报错烦的可以自己加个if语句修改下脚本。
  • 我们将K点增加到331,继续计算。

查看计算结果

查看结果的时候,有以下三个主要的方面:

1) 看能量

  • 可以使用VTST的nebresults.pl 的小脚本,也可以使用自己写的小脚本(ta.sh)。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ ls
00 01 02 03 04 05 06 07 08 09 exts.dat INCAR job_sub KPOINTS mep.eps movie movie.vasp neb.dat nebef.dat POTCAR slurm-1133307.out spline.dat vaspgr vasprun.xml
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ nebresults.pl

Unziping the OUTCARs ... done
Do nebbarrier.pl ; nebspline.pl
Do nebef.pl
Do nebmovie.pl
Do nebjmovie.pl
Do nebconverge.pl

Forces and Energy:
0 0.023782 -321.151800 0.000000
1 0.007650 -321.196900 -0.045100
2 0.009671 -321.160900 -0.009100
3 0.010446 -321.128300 0.023500
4 0.014498 -321.106600 0.045200
5 0.019054 -321.096100 0.055700
6 0.019424 -321.094300 0.057500
7 0.013325 -321.111700 0.040100
8 0.005225 -321.156900 -0.005100
9 0.011589 -321.098500 0.053300

Extremum 1 found at image 0.910082 with energy: -0.046551
Extremum 2 found at image 5.724723 with energy: 0.057913
Extremum 3 found at image 8.102362 with energy: -0.007602
Extremum 4 found at image 8.999187 with energy: 0.053367

iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ ta.sh
00 -321.15187450
01 -321.19698136
02 -321.16098546
03 -321.12832705
04 -321.10664059
05 -321.09613563
06 -321.09434359
07 -321.11177677
08 -321.15691232
09 -321.09850728

  • 看到能量的时候,前面我们粗算的时候,已经知道能量随着路径的变化是怎么样的,这里我们就脑子自动画图,从01到08,看着能量慢慢上升,在05,06左右的时候达到最高点,然后继续下降。

  • 注意:我们这是接着前面一节进行的计算, 00 和 09中的能量都是551K点下的能量,因此得到的能量随路径的变化是下面这样子的。如果看到这样的图片,不要慌张,跳过第一个和最后一个点,直接看中间的部分即可。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ grep irre */OUTCAR
    00/OUTCAR: Found 9 irreducible k-points:
    01/OUTCAR: Found 5 irreducible k-points:
    02/OUTCAR: Found 5 irreducible k-points:
    03/OUTCAR: Found 4 irreducible k-points:
    04/OUTCAR: Found 4 irreducible k-points:
    05/OUTCAR: Found 4 irreducible k-points:
    06/OUTCAR: Found 4 irreducible k-points:
    07/OUTCAR: Found 4 irreducible k-points:
    08/OUTCAR: Found 5 irreducible k-points:
    09/OUTCAR: Found 13 irreducible k-points:

2) 看计算有没有正常结束。

主要是查看:i) NEB跑了多少步,ii) 自己设置了多少步,iii) 通过OUTCAR中结构收敛的特定关键词。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ ls
00 02 04 06 08 exts.dat job_sub mep.eps movie neb.dat POTCAR spline.dat vasprun.xml
01 03 05 07 09 INCAR KPOINTS mep.png movie.vasp nebef.dat slurm-1133307.out vaspgr
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ tail 03/OSZICAR -n 5
RMM: 2 -0.318010871323E+03 -0.48063E-06 -0.15546E-06 134 0.370E-03 0.397E-03
RMM: 3 -0.318010871353E+03 -0.29988E-07 -0.12836E-08 21 0.273E-03 0.402E-03
RMM: 4 -0.318010871104E+03 0.24867E-06 0.00000E+00 17 0.260E-03 0.157E-03
RMM: 5 -0.318010871201E+03 -0.96799E-07 -0.11232E-08 18 0.266E-03
14 F= -.32114179E+03 E0= -.32112833E+03 d E =-.102536E-03
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ grep NSW INCAR
NSW = 200
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$ grep 'reached re' 03/OUTCAR
reached required accuracy - stopping structural energy minimisation
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex78$

当然你可以写一个小脚本如下来避免每次都重读在终端里面打出来。

1
2
3
4
5
#!/usr/bin/env bash 

tail 03/OSZICAR -n 5
grep NSW INCAR
grep 'reached re' 03/OUTCAR

也可以设置下alias来判断:

1
alias check='tail 03/OSZICAR -n 5 && grep NSW INCAR && grep "reached re" 03/OUTCAR '

有能力的可以写一些更高级点的,加些if,for语句智能判断NEB有没有算完,有没有收敛。

3) 查看NEB中各个IMAGE中的结构有没有跑乱,散架,有没有物理化学意义。前面我们讲过了几种办法,本节我们简单复习一下,就不再啰嗦了。

i) 通过vaspkit 结合 vmd实现看动画的效果(Windows用户)

ii)使用ASE 和 p4vasp 批量打开所有的IMAGES中的CONTCAR,挨个查看。

小节

本节,我们主要是在前面计算的基础上,把gamma点的输出结果备份了一下,在此基础上,增加K点(提高精度)继续优化。最近瞎忙了很多事情,让很多人等急了,给大家拜个晚年,祝大伙猪年都学会麻溜滴算过渡态。

通过job-ID快速进入计算目录,可以有效地避免狂输cd命令,提高我们的计算效率。本节我们通过在国科智算超算中心的演示,介绍一个Slurm任务管理系统,保存任务ID,计算目录,并快速根据ID进入计算目录的方法。

核心思想:

有很多种方法可以实现:通过Job-ID进入计算目录的办法。这里只介绍本节的思路。

  • 通过squeue 命令获取用户所有的任务的ID;
  • 通过scontrol 命令获取所有ID对应的计算目录
  • 将计算ID,以及对应的目录保存到一个txt文件中。
  • 通过在~/.bashrc中使用scource命令,激活bash脚本中进入目录的操作。

之所以将计算目录保存到txt文件中,是因为如果计算已经完成的话,则不可以通过scontrol命令来获取目录,这时候,我们就可以在txt文件中查找相应的信息。

准备流程

按照流程操作,每一步都很重要,漏掉的话,可能会导致功能无法实现。

1 下载脚本:

ent.shjob_check.py复制到~/bin 目录下,赋予可执行权限。

1
2
gkzshpc101@login02:~/bin$ chmod  u+x ent.sh
gkzshpc101@login02:~/bin$ chmod u+x job_check.py

2 在bin目录下创建一个:job-check的目录,以及在该目录中创建一个空的job_list.txt文件,用来存储我们的任务信息。

1
gkzshpc101@login02:~/bin$ mkdir job_check && touch job_check/job_list.txt

3 修改job_check.py脚本中

1)13行中的账户信息,gkzshpc101改成你自己的账户

2) 32行中的那个目录信息,改成下面目录输出的。

1
2
3
gkzshpc101@login02:~/bin$ cd job_check 
gkzshpc101@login02:~/bin/job_check$ pwd
/public1/home/gkzshpc101/bin/job_check
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#!/usr/bin/env python
'''
This script is used to
1) record the job_ID and directories
2) print the job path from the job_ID
'''
from subprocess import Popen, PIPE
import numpy as np
import sys

def get_id():
list_j = [] # list for the job_ID
process = Popen(['squeue', '-lu', 'gkzshpc101'], stdout=PIPE, stderr=PIPE)
stdout, stderr = process.communicate()
list_out = stdout.split('\n')[2:]
for i in range(0, len(list_out)-1):
list_j.append(list_out[i].split()[0])
return list_j

def get_dir(job_id):
job_dir = None
command = 'scontrol show job ' + job_id
process1 = Popen(command, shell = True, stdout=PIPE, stderr=PIPE)
stdout1, stderr1 = process1.communicate()
for i in stdout1.split('\n'):
if 'WorkDir' in i:
#job_dir.append(i.split('=')[1])
job_dir = i.split('=')[1]
return job_dir

id_dict = {}
with open('/public1/home/gkzshpc101/bin/job_check/job_list.txt', 'a+') as file_in:
lines = file_in.readlines()
for line in lines:
key = line.split(':')[0]
value = line.split(':')[1]
id_dict[key] = value

list_j = get_id()

for i in list_j:
if i not in id_dict.keys():
id_dict[i] = get_dir(i)
file_in.write(i + ':' + get_dir(i) + '\n')

job_id = sys.argv[1]

for i in id_dict.keys():
if job_id in i:
print(id_dict.get(i).strip())

1
2
3
4
#!/usr/bin/env bash
pwd_work=$(job_check.py $1)
echo $pwd_work
cd $pwd_work

4 修改~/.bashrc文件,添加下面这一行:

1
alias ent='source ~/bin/ent.sh '

source 一下~/.bashrc文件: . ~/.bashrc

运行实例:

1
2
3
4
5
6
7
8
9
10
11
12
13
gkzshpc101@login02:~/test_ncore/test1/14$ ls
INCAR KPOINTS POSCAR POTCAR vasp.sh
gkzshpc101@login02:~/test_ncore/test1/14$ sbatch vasp.sh
Submitted batch job 14999
gkzshpc101@login02:~/test_ncore/test1/14$ squeue
JOBID PARTITION NAME USER ST TIME NODES NODELIST(REASON)
14999 operation vasp gkzshpc1 R 0:02 1 c0077
gkzshpc101@login02:~/test_ncore/test1/14$ cd
gkzshpc101@login02:~$ ent 999
/public1/home/gkzshpc101/test_ncore/test1/14
gkzshpc101@login02:~/test_ncore/test1/14$ ls
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT INCAR KPOINTS OSZICAR OUTCAR PCDAT POSCAR POTCAR REPORT vasp.out vasprun.xml vasp.sh WAVECAR XDATCAR
gkzshpc101@login02:~/test_ncore/test1/14$

国科智算提交VASP任务

今天QQ群里,有人问怎么在国科智算的超算中心提交VASP的任务。本着一言不合就写教程的态度,这一节我们就看下VASP的任务是怎么提交的。提交任务的脚本见群文件:vasp_qiangli.sh。下载后重命名成vasp.sh即可。想试用的,购买机时的可以加QQ群:608307988咨询一下。

Slurm

国科智算上系统采用的是slurm 任务调度工具。Simple Linux Utility for Resource Management,它是一个用于LinuxUnix 内核系统的免费、开源的任务调度工具,被世界范围内的超级计算机和计算机群广泛采用。它提供了三个关键功能。第一,为用户分配一定时间的专享或非专享的资源(计算机节点),以供用户执行工作。第二,它提供了一个框架,用于启动、执行、监测在节点上运行着的任务(通常是并行的任务,例如MPI,第三,为任务队列合理地分配资源。 大约60%的500强超级计算机上都运行着Slurm,包括2016年前世界上最快的计算机天河-2。 以上是来自维基百科的解释,具体的大家可以浏览Slurm 的官网:Slurm Workload Manager

Sbatch

在这个调度系统中,提交任务时我们需要用到命令: sbatch : https://slurm.schedmd.com/sbatch.html

sbatch 命令后面要跟一堆的参数,比如计算时间,节点数,邮箱,队列,调用的环境变量,任务名称等等。但这些信息通过命令直接输入又有些麻烦,所以我们把它们放到一个脚本里面,免得每次都重新输入一大长串的内容。而这个脚本,也就是我们本节内容的主角: vasp.sh

首先,我们浏览sbatch的详细参数: https://slurm.schedmd.com/sbatch.html

然后根据这些参数,我们就可以创建一个vasp.sh脚本了。

vasp.sh

下面就是vasp.sh 的主要内容:有2种写法,这两种写法也可以混着用,不影响。主要还是要参考sbatch的使用说明,需要什么就按照格式填写相应的内容。

写法(一)

1
2
3
4
5
6
7
8
9
10
11
12
#!/bin/bash
#SBATCH -J BigBro ## Job Name
#SBATCH -o %j.out ## standard output
#SBATCH -e %j.err ## standard error
#SBATCH -p operation ## Partition
#SBATCH -N 1 ## Number of nodes
#SBATCH --ntasks-per-node=28 ## Each node has 28 tasks
#SBATCH -t 02-23:57:25 ## time for your job: 2 d,23 h ,57 min and 23 s

module load mpi/intelmpi/2017.4.239
mpirun /public/software/apps/vasp/544/vasp.5.4.4/bin/vasp_std

写法(二)

1
2
3
4
5
6
7
8
9
10
11
#!/bin/bash
#SBATCH --job-name=BigBro ## Job Name
#SBATCH --output=%j.out ## standard output
#SBATCH --error=%j.err ## standard error
#SBATCH --partition=operation ## Partition
#SBATCH --nodes=1 ## Number of nodes
#SBATCH --ntasks-per-node=28 ## Each node has 28 tasks
#SBATCH --time=02-23:57:25 ## time for your job: 2 d,23 h ,57 min and 23 s

module load mpi/intelmpi/2017.4.239
mpirun /public/software/apps/vasp/544/vasp.5.4.4/bin/vasp_std

上面每一行的含义,大师兄都注释出来了。

从第二行往下一次为:

  • 任务的名字,
  • VASP的标准输出,
  • 错误输出,
  • 任务运行的分区,
  • 使用的节点,
  • 每个节点的的tasks数目。这个tasks可以理解为每个节点的核数,国科智算的新机器每个节点是28个核。
  • 以及你的计算所用的时间。

最后两行为:我们调用的环境变量以及运行vasp程序。

需要注意的是,脚本里面的内容比如-N(修改任务所需的节点数目)、-J(修改任务的名字)这些我们频繁更换的,可以从脚本里面拿出来,在命令中运行。

1
sbtach -N 2 -J test vasp.sh 

本人自己的计算,一般28个核就够,任务名字也懒得修改。就把-N,-J写到vasp.sh中了。大家根据自己的任务特点自动修改就行了。下面我们具体演示一下。

实例操作1:

我们要在operation分区,运行一个VASP的单点计算,任务名称为:single, 使用2个节点,限制时间为2个小时;那么脚本的修改以及任务提交如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
gkzshpc101@login02:~/ex-A11/test_single$ 
gkzshpc101@login02:~/ex-A11/test_single$ ls
INCAR KPOINTS POSCAR POTCAR vasp.sh
gkzshpc101@login02:~/ex-A11/test_single$ cat vasp.sh
#!/bin/bash
#SBATCH -J single
#SBATCH -o out.%j
#SBATCH -e err.%j
#SBATCH -p operation
#SBATCH -N 2
#SBATCH --ntasks-per-node=28
#SBATCH -t 02:00:25

module load mpi/intelmpi/2017.4.239
mpirun /public/software/apps/vasp/544/vasp.5.4.4/bin/vasp_std
gkzshpc101@login02:~/ex-A11/test_single$ sbatch vasp.sh
Submitted batch job 35525
gkzshpc101@login02:~/ex-A11/test_single$ squeue
JOBID PARTITION NAME USER ST TIME NODES NODELIST(REASON)
35525 operation single gkzshpc1 R 0:09 2 c[0032-0033]
gkzshpc101@login02:~/ex-A11/test_single$ ls
CHG CHGCAR CONTCAR DOSCAR EIGENVAL err.35525 IBZKPT INCAR KPOINTS OSZICAR out.35525 OUTCAR PCDAT POSCAR POTCAR REPORT vasprun.xml vasp.sh WAVECAR XDATCAR
gkzshpc101@login02:~/ex-A11/test_single$

实例操作2:

我们在Operation分区,运行一个VASPCI-NEB的过渡态计算任务,任务名称为:NEB,使用4个节点,限制时间为12个小时;那么脚本的修改以及任务提交如下:插了7个点,用112个核算,16个核算一个点。下面有点长,文字描述就到此结束,自己慢慢看,希望对大家有所帮助。注意:脚本里面,我们换成编译了VTSTvasp 5.4.1版本。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
gkzshpc101@login02:~/ex-A11/test_neb$ ls
00 01 02 03 04 05 06 07 08 FS INCAR IS KPOINTS POTCAR vasp.sh
gkzshpc101@login02:~/ex-A11/test_neb$ cat vasp.sh
#!/bin/bash
#SBATCH -J NEB
#SBATCH -o out.%j
#SBATCH -e err.%j
#SBATCH -p operation
#SBATCH -N 4
#SBATCH --ntasks-per-node=28
#SBATCH -t 12:00:25

module load mpi/intelmpi/2017.4.239
mpirun /public/software/apps/vasp/541_neb/vasp.5.4.1/bin/vasp_std
gkzshpc101@login02:~/ex-A11/test_neb$
gkzshpc101@login02:~/ex-A11/test_neb$ sbatch vasp.sh
Submitted batch job 35526
gkzshpc101@login02:~/ex-A11/test_neb$ squeue
JOBID PARTITION NAME USER ST TIME NODES NODELIST(REASON)
35526 operation NEB gkzshpc1 R 0:07 4 c[0056-0057,0060-0061]
gkzshpc101@login02:~/ex-A11/test_neb$ ls
00 01 02 03 04 05 06 07 08 err.35526 FS INCAR IS KPOINTS out.35526 POTCAR vasprun.xml vasp.sh
gkzshpc101@login02:~/ex-A11/test_neb$ ls *
err.35526 FS INCAR IS KPOINTS out.35526 POTCAR vasprun.xml vasp.sh

00:
POSCAR

01:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT WAVECAR XDATCAR

02:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

03:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

04:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

05:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

06:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

07:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR

08:
CHG CHGCAR CONTCAR DOSCAR EIGENVAL IBZKPT OSZICAR OUTCAR PCDAT POSCAR REPORT stdout WAVECAR XDATCAR
gkzshpc101@login02:~/ex-A11/test_neb$
gkzshpc101@login02:~/ex-A11/test_neb$ head -n 10 out.35526
running on 112 total cores
each image running on 16 cores
distrk: each k-point on 16 cores, 1 groups
distr: one band on 1 cores, 16 groups
vasp.5.4.1 05Feb16 (build May 24 2018 19:36:47) complex

前面的准备工作完成了,INCARKPOINTSPOSCAR(IMAGES), POTCAR也检查完了。剩下的就是准备脚本提交文件了。本节主要是在天河II号超算中心上给大家简单示范一下:

  • 1)准备脚本,提交任务;
  • 2)过渡态任务运行时候的查看;
  • 3) VTST脚本nebresults.pl的安装和使用。

提交任务

首先通过命令操作,熟悉下天河II号超算中心提交NEB计算的流程。

1
2
3
4
5
6
7
8
9
10
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ ls
00 01 02 03 04 05 06 07 08 09 INCAR job_sub KPOINTS POTCAR
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ cat job_sub
#!/bin/bash
export LD_LIBRARY_PATH=/THFS/opt/intel/composer_xe_2013_sp1.3.174/mkl/lib/intel64:$LD_LIBRARY_PATH
yhrun -p gsc -n 24 /THFS/opt/vasp/5.4.4_neb/vasp.5.4.4/bin/vasp_std
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77/$ yhbatch -p gsc -N 2 -J test job_sub
Submitted batch job 1004311
  • 每个节点有24核,我们使用了2个节点,共48核,来计算这个NEB任务。
  • 共有8个IMAGES,也就是6个核算一个IMAGE。
  • 我们调用的VASP版本是:/THFS/opt/vasp/5.4.4_neb/vasp.5.4.4/bin/vasp_std
  • 本例子下载链接:https://pan.baidu.com/s/1shjrJsRuDNYGEd2qDqvdFQ 提取码:k3pp

查看NEB计算

NEB一般收敛的很慢,需要跑很多离子步才会收敛,所以大家要心里有个准备,这也是为什么我们先用gamma点粗算一下,再用高密度的K点计算的原因,总之就是节约时间。虽然NEB要花很长时间才收敛,但提交任务后,我们不能守株待兔似的等着NEB算完,而是要勤检查结果,因为NEB的计算中,也会经常出现结构跑乱的情况。那么该怎么检查呢? 一看能量,二看结构,三看能量和结构。

一看能量:

主要是因为在terminal下面,相对于打开可视化软件查看结构来说,我们通过命令提取每个IMAGEOUTCAR能量信息的操作更加方便些而已。更重要的其实还是结构的变化,也就是你的反应路径。查看能量有2个办法:

自己用脑子想象或者写脚本:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ ls
00 01 02 03 04 05 06 07 08 09 INCAR job_sub KPOINTS NEB.pdb POTCAR slurm-995085.out vasprun.xml
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ ta.sh
01 -298.00705385
02 -297.96462657
03 -297.92754864
04 -297.90142637
05 -297.88475532
06 -297.87681445
07 -297.89025285
08 -297.93374021
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ cat ~/bin/ta.sh
#!/usr/bin/env bash
for i in *; do
if [ -e $i/OUTCAR ]; then
echo -e $i "\t" $(grep ' without' $i/OUTCAR |tail -n 1 | awk '{print $7}')
fi
done
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$

用脑子想象:

ts.sh 脚本的输出中,从01的能量一直往下看,自己脑子里面自带一个xy的坐标系,将这些能量填到坐标系中,直至08。你会想象得到这样的一个曲线: 能量从01的-298.007慢慢上升到06的-297.87,然后再下降到08的-297.934。这个曲线的顶点再06处,也就是我们粗算得到的过渡态。

当然,你也写个小脚本提取每个IMAGE中的能量信息,然后画图即可。

使用VTST的脚本:nebresults.pl

前面Ex72中,我们讲过了怎么准备VTST的那些脚本,这一节,我们讲另外一个办法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
iciq-lq@ln3:/THFS/home/iciq-lq/bin$ cd vtst/
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst$ ls
vtstscripts.tgz
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst$ tar -zxvf vtstscripts.tgz
vtstscripts-937/
vtstscripts-937/pos2xyz.py
vtstscripts-937/sum_dos_np
vtstscripts-937/chg2cube.pl
vtstscripts-937/chgsplit.sh
vtstscripts-937/akmcprocess.pl
*
*
*
vtstscripts-937/nebbarrier.pl
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst$ ls
vtstscripts-937 vtstscripts.tgz
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst$ cd vtstscripts-937/
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst/vtstscripts-937$ ls
2con.py center.py diffcon.pl dymanalyze.pl dymseldsp.pl insplot.pl ....
  • 使用vim打开nebresults.pl文件,将58-71行注释掉,或者可以使用sed命令:
1
sed -i '58,71s/^/#/g' nebresults.pl

效果如下图:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
57 # Zip the OUTCARs again
58 #print 'Zipping the OUTCARs again ... ' ;
59 #$zip = $ENV{'VTST_ZIP'};
60 #if($zip eq ''){ $zip = 'gzip'; }
61 #
62 #$i = 0;
63 #$string = "00";
64 #while(chdir $string) {
65 # system "$zip OUTCAR";
66 # $i++;
67 # if($i < 10) { $string = "0$i"; }
68 # elsif($i < 100) { $string = "$i"; }
69 # chdir $dir;
70 #}
71 #print "done\n";
  • ~/.bashrc文件中设置脚本的目录:
1
2
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst/vtstscripts-937$ pwd
/THFS/home/iciq-lq/bin/vtst/vtstscripts-937
  • 复制上面pwd命令的目录,打开~/.bashrc 文件,添加这一行:
1
export PATH=$PATH:/THFS/home/iciq-lq/bin/vtst/vtstscripts-937
  • source一下~/.bashrc文件:
1
iciq-lq@ln3:/THFS/home/iciq-lq/bin/vtst/vtstscripts-937$ . ~/.bashrc
  • 进入计算的目录下,运行nebresults.pl命令,操作如下图:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ ls
00 01 02 03 04 05 06 07 08 09 INCAR job_sub KPOINTS NEB.pdb POTCAR slurm-995085.out vasprun.xml
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ nebresults.pl

No OUTCAR in 00
Unziping the OUTCARs ...
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ cp 01/OUTCAR 00
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ nebresults.pl

No OUTCAR in 09
Unziping the OUTCARs ...
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ cp 08/OUTCAR 09
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ nebresults.pl

Unziping the OUTCARs ... done
Do nebbarrier.pl ; nebspline.pl
Do nebef.pl
Do nebmovie.pl
Do nebjmovie.pl
Do nebconverge.pl

Forces and Energy:
0 0.003505 -298.009300 0.000000
1 0.003505 -298.009300 0.000000
2 0.008701 -297.969500 0.039800
3 0.015508 -297.932900 0.076400
4 0.018300 -297.905700 0.103600
5 0.016765 -297.887500 0.121800
6 0.011627 -297.877000 0.132300
7 0.007293 -297.885300 0.124000
8 0.003607 -297.929100 0.080200
9 0.003607 -297.929100 0.080200

Extremum 1 found at image 0.210729 with energy: 0.002874
Extremum 2 found at image 0.788081 with energy: -0.002910
Extremum 3 found at image 6.265931 with energy: 0.132792
Extremum 4 found at image 8.211450 with energy: 0.074781
Extremum 5 found at image 8.788801 with energy: 0.085561

iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$

上面的那个图在Ubuntu的终端里面,通过gs mep.eps命令可以直接打开。

解释:

1) 使用nebresults.pl脚本的时候,我们需要在0009的文件夹中分别放上初始结构和末态结构所对应的OUTCAR。上面操作中,大师兄把0108中的OUTCAR分别复制到了0009中。

为什么这么做呢?

答:偷懒。原因有2个:

A)之前优化结构都是用到的高密度的K点,直接把OUTCAR拿过来用的话,会出现能量差别很大。做出图来能量很奇怪的样子。如下面的例子所示:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
  iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77/high_k$ ta.sh
00 -321.15187450
01 -298.00935617
02 -297.96955082
03 -297.93296493
04 -297.90579729
05 -297.88754704
06 -297.87700139
07 -297.88531984
08 -297.92917858
09 -321.09850728
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77/high_k$ nebresults.pl

Unziping the OUTCARs ... done
Do nebbarrier.pl ; nebspline.pl
Do nebef.pl
Do nebmovie.pl
Do nebjmovie.pl
Do nebconverge.pl

Forces and Energy:
0 0.023782 -321.151800 0.000000
1 0.003505 -298.009300 23.142500
2 0.008701 -297.969500 23.182300
3 0.015508 -297.932900 23.218900
4 0.018300 -297.905700 23.246100
5 0.016765 -297.887500 23.264300
6 0.011627 -297.877000 23.274800
7 0.007293 -297.885300 23.266500
8 0.003607 -297.929100 23.222700
9 0.011589 -321.098500 0.053300

Extremum 1 found at image 0.000009 with energy: -0.000000
Extremum 2 found at image 6.265997 with energy: 23.275310

iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$

看到了没有,上面例子中0108 结构的能量很大。这是因为参考的能量是00的。对于同一个体系,K点大的话,绝对能量更负一些。能量差别太大,导致过渡态的能量变化都可以忽略掉了。

B)当然我们也可以分别对0009对应对的机构算个单点。生成各自对应的OUTCAR,然后再使用nebresults.pl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ ta.sh
00 -298.03233582
01 -298.00705385
02 -297.96462657
03 -297.92754864
04 -297.90142637
05 -297.88475532
06 -297.87681445
07 -297.89025285
08 -297.93374021
09 -297.97265721
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$
iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$ nebresults.pl

Unziping the OUTCARs ... done
Do nebbarrier.pl ; nebspline.pl
Do nebef.pl
Do nebmovie.pl
Do nebjmovie.pl
Do nebconverge.pl

Forces and Energy:
0 0.100196 -298.032300 0.000000
1 0.003505 -298.009300 0.023000
2 0.008701 -297.969500 0.062800
3 0.015508 -297.932900 0.099400
4 0.018300 -297.905700 0.126600
5 0.016765 -297.887500 0.144800
6 0.011627 -297.877000 0.155300
7 0.007293 -297.885300 0.147000
8 0.003607 -297.929100 0.103200
9 0.084593 -297.972600 0.059700

Extremum 1 found at image 0.057295 with energy: -0.000155
Extremum 2 found at image 6.265997 with energy: 0.155771

iciq-lq@ln3:/THFS/home/iciq-lq/LVASPTHW/ex77$

仔细观察,你会发现,00010809 的能量差别基本不大。所以上面我们的偷懒做法也是可取的。

敲黑板:

一般在粗算的时候,不用考虑那么多细节,怎么粗怎么弄。所以,我们前面讲的这个懒办法是在粗算情况下的操作。由于高K点的计算我们在优化初末态结构的时候已经有了,所以提高精度的时候就可以直接用了。

本节要求

1) 会提交NEB的计算任务

2) 会使用自己写的或者nebresults.pl脚本查看能量信息。

下一节,我们学习怎么分析这些能量。