VASP计算金刚石硬度(体积弹性模量)
课题的意义就不说了，bay拿这个体系练练手学习一下VASP感觉结果不错，网上也鲜有详细的教程，这里写一点
一是给象bay这种接触VASP不到一周的人介绍一下用法
二是请教熟悉者分析和优化一下模拟参数
金刚石的立方体cell如图所示，每条边长l=3.567Ang，一个cell内有8个原子。
其primary cell可用三条夹角60度的基失表示:(a,a,0),(a,0,a),(0,a,a),
其中a＝l/2，内部包含2个原子，分数坐标为(0,0,0),(1/4,1/4,1/4).
text的图示都没对齐，这里传一个文本。
本文内容遵从CC版权协议 转载请注明，更新以本网页为准。
1 这里首先让VASP跑起来
运行VASP需要在所在目录内准备4个文件，
POSCAR：体系结构，包含cell大小和原子坐标
INCAR：绝大部分模拟参数
KPOINTS：K Mesh
POTCAR：赝势文件，这个别人作好了我们可以直接拿来用，不需要自己写
POSCAR：对着上面第二段文字可以大概猜出Diamond
1.7835
0.01.01.0
1.00.01.0
1.01.00.0
2
Direct
0.000.000.00
0.250.250.25复制代码INCAR：System = diamond, Opt cell
ENCUT = 400
GGA = 91
ISMEAR = -5
PREC = Accurtae
IBRION = 2
ISIF = 7
NSW = 50复制代码ENCUT = 400:截断能为400eV
GGA = 91 :对赝势进行GGA计算，采用PW91泛函，比较常用的方法
PREC = Accurtae:较高的精度
IBRION;ISIF;NSW:对cell进行优化
KPOINTS：使用8*8*8的K点数目K MESH
0
Monkhorst
8 8 8
0 0 0复制代码POTCAR 一般是每个元素一个文件，注意泛函和INCAR中指定的相一致
进入存放这四个文件的目录，允许vasp，就可以进行计算。体系较小，数秒即可完成，不需要并行。
完成后会生存很多文件，这里我们只要
OUTCAR：详细的输出文件
提取最终的能量值$ grep TOTEN OUTCAR | tail -n 1
freeenergy TOTEN= -18.201948 eV复制代码和CONTCAR(优化后的坐标文件)，前半部分的格式同输入文件POSCAR相同Diamond
1.78350000000000
0.00000000000000000.99904263112518790.9990426311251879
0.99904263112518790.00000000000000000.9990426311251879
0.99904263112518790.99904263112518790.0000000000000000
Direct
0.00000000000000000.00000000000000000.0000000000000000
0.25000000000000000.25000000000000000.2500000000000000
...复制代码1.7835*0.999即为优化后的cell尺寸(平衡点)。INCAR中的ISIF=7表示仅优化cell不优化原子位置，后面两行未变 2 选用不同的参数，对模拟条件进行优化。首先是ENCUT的大小
计算在平衡点附近，不同cell大小的能量值，观测起连续性
删除INCAR中的IBRION;ISIF;NSW三行，连续改变cell尺寸，从OUTCAR文件提取能量值
为简化操作，这里写了一个脚本run.sh#!/bin/bash
for i in 1.781 1.782 1.783 1.784 1.785 1.786 1.787 1.788 1.789 1.79 1.791 1.792 1.793 1.794
do
cat > POSCAR > /dev/null
E=`grep TOTEN OUTCAR | tail -n 1 | awk '{print($5)}'`
echo $i $E
done复制代码运行方法：
./run.sh | tee EC400.out
cell尺寸和能量值都会写到EC400文件中，每行一个。修改INCAR中的ENCUT为300和500，重复操作。
结果如下，用gnuplot作的图，至少我用起来很舒服。
更多命令行下的作图方法参考之前的这个帖子 http://www.mdbbs.org/thread-24037-1-1.htmlgnuplot> plot "EC500"
-18.1924 ---------------------------------------------------------
A "EC500K8" AA
-18.1926
| |
-18.1928
| |
-18.193 A A
| |
-18.1932
| |
-18.1934 A
| A |
-18.1936
| |
-18.1938 A A
| |
-18.194
| A A |
-18.1942 A A
AA
-18.1944 ---------------------------------------------------------
1.781.782 1.7841.786 1.788 1.79 1.7921.794
gnuplot> plot "EC400"
-18.2 ---------------------------------------------------------
"EC400K8" A
| A
-18.2005
| |
|A A|
|A |
-18.201
| |
| |
-18.2015 A A
| A A |
| |
-18.202
| |
| A |
| A |
-18.2025 A AAA
| |
-18.203 ---------------------------------------------------------
1.781.782 1.7841.786 1.788 1.79 1.7921.794
gnuplot> plot "EC300"
-18.334 ----------------------------------------------------------
"EC300K8" A
-18.335 A
| |
| |
-18.336
| |
-18.337 A
| A A |
| AAA |
-18.338
| |
-18.339 A A
| |
| A |
-18.34
|A A |
-18.341 A
| |
A
-18.342 ----------------------------------------------------------
1.78 1.7821.784 1.7861.7881.791.792 1.794复制代码可以看到ENCUT＝300时候，结果已经乱了。400可用，条件允许时候用500吧。 K Mesh的选取
这里我们正式计算一个体积弹性模量bulk modulus，用精度较高的ENCUT＝500和K mesh 8*8*8
那个run.sh脚本改成这个样子#!/bin/bash
for i in 1.71 1.72 1.73 1.74 1.75 1.76 1.77 1.78 1.79 1.8 1.81 1.82 1.83 1.84
do
cat > POSCAR > /dev/null
E=`grep TOTEN OUTCAR | tail -n 1 | awk '{print($5)}'`
echo $i $E
done复制代码运行 run.sh | tee run.out
可以得到每个cell尺寸对应的体系能量值，根据Murnaghan状态方程 http://en.wikipedia.org/wiki/Bir ... n_equation_of_state
E=E0(B/2V)(V-V0)^2 (*)
即可拟合出B的数值，即为体系的Bulk modulus，体积随压强的变化性能。其中V0，E0是平衡点的体积和压强值。
该系列点近似是一个二次曲线，像下面这个图一样。先求出这个曲线的对称轴和最低点吧，也就是V0和E0，图里的B点是bay自己求出的，方法见下-17.9 ---A-------------------------------------------------------
"run.out" A
| |
-17.95
| |
| A |
| |
-18
| |
| A |
-18.05
| |
| A
-18.1 A
| |
| A A|
| |
-18.15 A A
| A |
AA A
-18.2 -------------------------------------BA--------------------
1.7 1.72 1.741.761.781.81.821.84
l0= 1.78755 Ang, e0= -18.194342 eV复制代码有个现成的murn.f代码 http://www.fhi-berlin.mpg.de/th/fhi98md/Murn/readme_murn.html
可以根据上面的run.out的数值自动拟合出B。
一个完整的输入文件如下，命名为inp.m2
0.25
1.7 1.9 5
14
1.71 -17.916690
1.72 -17.988714
1.73 -18.050681
1.74 -18.102208
1.75 -18.141745
1.76 -18.171599
1.77 -18.191770
1.78 -18.201224
1.79 -18.202421
1.8 -18.193810
1.81 -18.180601
1.82 -18.156983
1.83 -18.127412
1.84 -18.090229复制代码第一行为能量单位，2表示eV
第二行为primarycell同cell的体积比，这里为0.25
第三行为拟合的范围，从1.7到1.9取50个点
第四行为14个数据，接下来是这14个cell尺寸vs能量
运行方法
murn
这个程序号称可以直接计算出B，好像有些问题，这里只拿它拟合V0。
也就是在输出文件中找这一行
alat= 1.78755 b0= ...........
那个1.78755就是V0，上面那个图的对称轴
拿这个数值带入POSCAR，计算得到这时的TOTEN，就是图中的最低点，Murnaghan方程中的E0。 将Murnaghan方程E=E0(B/2V)(V-V0)^2 变换成
B＝2V0*(E-E0)/(V-V0)^2 (**)
逐点带入即可求出各个B
一个简单的脚本如下
awk '{print($1*$1*$1*8,$2*4)}' run.out > VE
新生成的VE第一列为立方体cell的体积，第二列为8个原子的能量值。40.0017 -71.6668
40.7076 -71.9549
41.4217 -72.2027
42.1442 -72.4088
42.875 -72.567
...复制代码v=45.6225198
e=72.811844
eval "awk '{print(2*\$v*(\$2$e)/(\$1-$v)/(\$1-$v)*160.2)}' VE" > B复制代码变量v和e是平衡点的V0和E0，只是把e的负号去掉了.{}里面就是方程(**)，160.2是eV/Ang^3到GPa的单位换算.
这时会得到一系列的B,bay作出的结果是463.908,459.601,455.241,450.847,446.447,441.918,436.713,429.111,485.134,431.275,425.426,420.599,416.31,412.11
490 -------------------------------------------------------------
A "B0" A
480
| |
| |
470
| |
460 A
|A |
| A |
450 A
| A |
440 A
| A |
| |
430 A A
| A |
420 A
| A |
A
410 -------------------------------------------------------------
0 2 4 6 8 101214复制代码两端的数值会因偏离Murnaghan 方程而带来误差，中间的会因大数相减带来误差，总的来说各个点都还规矩
可以拟合直线，也可以在平衡点附近求平均。这几个点的平均值约为441 GPa，同实验值442 GPa极为接近。
Linux下没有找到特别好用的轻量级拟合软件，win下可以用用这个
结果汇总
ENCUTK :B
500 8*8*8: 441.0
400 8*8*8: 435.7
500 4*4*4: 437.7
400 4*4*4: 438.2
这个方法，即使在不高的水平下，结果也是比较精确的。 比较详细！ 很详细，对初学者帮助很大。 :handshake
页:
[1]