在材料力学强度的表征测试中,理想拉伸曲线是其中典型的一种方式。从拉伸曲线中我们可以获得弹性模量,屈服强度等等有用的力学系数。使用第一性原理模拟理想拉伸(不局限于拉伸,但要注意,不要将宏观应力应变与DFT计算的理想应力应变混淆)曲线当然是一种很fancy的方法,但是目前我没有看到一个自动化实现这个过程的脚本。我用bash写了一个idealdeform.sh脚本,可以自动化实现理想拉伸或者剪切过程,目前主要用于VASP。不要问我为什么用bash,理由详见上篇博文(vasp定压计算脚本vaspeqstress.sh使用教程)。
![e2e4521681e08ac4941c189e7d2c59da.png](https://i-blog.csdnimg.cn/blog_migrate/6ef83dd6ef73db4ba930297d3dfef32d.jpeg)
Nanoscale, 2017, 9, 850–855
还是先讲讲原理吧。脚本原理相当简单,我主要先讲一下两类应变和应力。我们知道,我们平时主要用的是工程应变和工程应力。这类应变应力定义如下:
![dbfaaea523324ba23b6d0c49ae5e7a42.png](https://i-blog.csdnimg.cn/blog_migrate/c21cc0d5c6d0c1a27cb1d374b9206b8e.jpeg)
![db0e01eecc9dca4a1534e0986974924c.png](https://i-blog.csdnimg.cn/blog_migrate/7485b7bf05b26f892aae753904071401.jpeg)
各项符号我就不过多解读了,这个应该是基础的东西。我们注意到,l_0 和A_0都是样品在未拉伸时的初始状态,这也就导致工程应变和工程应力并非材料实际所受的状态,因此我们日常看到的实际样品的工程应力应变曲线往往都是下面这种低碳钢形状的。
![e0247aae4943e7357c8946f8095029c2.png](https://i-blog.csdnimg.cn/blog_migrate/0e93a26f1f3d0ece8d25ce8fbf839527.jpeg)
https://www.scienceabc.com/innovation/what-is-the-stress-strain-curve.html
材料在B点开始屈服,D点开始发生颈缩,E点发生断裂。但实际上材料发生颈缩后,材料继续加工硬化,真实应力实际上是继续上升,直到断裂。为了更加准确地描述这个过程,有时候大家也会使用真实应变和真实应力来绘制曲线。下面我们来简单推导一下两者地表达式。
考虑塑性变形的不可压缩性,我们有:
![6d0645b6b7a44c141a082283940f4bf3.png](https://i-blog.csdnimg.cn/blog_migrate/b3d88fe1a803a45798538e38282436d5.jpeg)
我们将当前面积A的定义代入到真实应力的定义式中即得真实应力和工程应力的换算关系。
![4247ab0e4737d316e4601906f2c3aebf.png](https://i-blog.csdnimg.cn/blog_migrate/00ef4f3ee9d1a005f8da8f17e12cd217.jpeg)
真实应变就更加简单了,我们只需要对每个瞬时应变进行积分即可
![ccb2738790cbb100a293d6224d5886dd.png](https://i-blog.csdnimg.cn/blog_migrate/e229dc01bf38ec238555398e674ea1e0.jpeg)
但实际上这种真实应变与工程应变的换算关系仅存在于颈缩初期之前,因为后期的时候材料的不可压缩的假设已经崩溃。因此实际材料变形曲线还需进一步校正,不过本文就不再过多讨论了。再次强调,我只是介绍了两类应力应变,请不要将宏观应力应变和理想应力应变联系起来。
下面我开始讲一下脚本的使用方法,原理就不细讲了,根本在于假设应变过程足够慢,可以看作准静态应变,这样我们只需按一定应变间隔施加应变矩阵,优化晶胞和原子即可。首先我们需要准备优化晶胞所需要的VASP四个输入文件。注意,POSCAR的原子坐标必须是分数坐标 。INCAR中ISIF可以用4,但请在OPTCELL中关闭你变形方向的弛豫 (OPTCELL是什么?自己去DET群文件搜搜吧)。最好不要用3,晶格的体积发生变化会导致曲线失真。然后将idealdeform.sh放到同一文件夹。我测试的体系是fcc Ni, a为3.5058Å。
现在我们来修改脚本中相关参数。
![098e1c6733d80003a3f5b19564d5b6b8.png](https://i-blog.csdnimg.cn/blog_migrate/cfe15c18e71be22083e04039c0c24c35.jpeg)
orientation表示施加应变的方式,你可以选择XX, YY, ZZ, XY, XZ, YZ六种变形方式中的一种,前三种为拉伸,后三种为剪切。initial为施加的初始应变,step表示每一步应变的间隔,num表示从初始应变开始按照step间隔一共连续施加了多少次应变。图中的例子就表示对fcc Ni沿着X轴拉伸,初始应变为0,每隔0.01施加应变进行一次拉伸,一共施加了100次应变,最后施加的应变为1.000(100%)。mpiexec表示你的运行语句,由你的系统决定。最后直接运行脚本即可。
![6afc8907109b92b9c899c9a2bb967717.png](https://i-blog.csdnimg.cn/blog_migrate/5a6f325ba999c3fdf38923b172702684.jpeg)
脚本会告诉你目前进行的进度,我们只需静待完成。计算完成后,程序默认会用gnuplot画一个工程应力应变曲线。
![c05e94f4b05e9d8fc50171fbc049b068.png](https://i-blog.csdnimg.cn/blog_migrate/aa86280dc36f9a41016491a3a05c013c.jpeg)
单位为GPa。我已经重新修改过单位,正号就表示拉应力。如果你觉得这样子看很辣眼睛,你可以将plot.sh也复制到当前文件夹,然后运行。
![bb9551010ecd97960577a23c0cd6b7f0.png](https://i-blog.csdnimg.cn/blog_migrate/c824ebc3afd869698404db7284db36fa.jpeg)
脚本运行完后,工程应力应变曲线的数据会被保存到engineeringstressstrain.all中,真实应力应变曲线的数据会被保存到truestressstrain.all中。如果某一步中出现无法收敛的情况,程序会停止。你应该在OUTCAR中查找原因,并从当前阶段继续运算。
下面说几个注意事项,非常重要:
大部分时间你并不需要用到真实应变,尤其当你的体系存在表面的时候,你不用去管真实应力应变曲线
DFT获得的应力应变曲线跟宏观应力应变曲线没有任何可比性,使用者应该清楚自己使用这个脚本获得的应力应变曲线到底应该如何理解。宏观应力应变来源于位错运动,跟DFT结果完全不是一个概念,我之所以在开头描述了一堆主要是为了介绍真实应力应变。
上述脚本在使用中如果有任何问题请联系18709821294@outlook.com,请大家关注ponychen的GitHub主页:
https://github.com/ponychen123/Vasptools
上面脚本和教程下载链接:
http://lanzous.com/u/dft_family
MS切面及构建根号表面经验
MATLAB做晶体结构图(固体物理)
固体物理学习笔记-能带理论(一)
量子力学笔记—什么是狄拉克符号?