课题的意义就不说了，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.0  1.0  1.0
1.0  0.0  1.0
1.0  1.0  0.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.0000000000000000  0.9990426311251879  0.9990426311251879
0.9990426311251879  0.0000000000000000  0.9990426311251879
0.9990426311251879  0.9990426311251879  0.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" A  A
-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        
             A  A      
-18.1944 ---------------------------------------------------------
   1.78  1.782 1.784  1.786 1.788 1.79 1.792  1.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 A  A  A        
    |                     |
                     
-18.203 ---------------------------------------------------------
   1.78  1.782 1.784  1.786 1.788 1.79 1.792  1.794
gnuplot> plot "EC300"
-18.334 ----------------------------------------------------------
                 "EC300K8" A  
-18.335                       A
   |                      |
   |                      |
-18.336                       
   |                      |
-18.337                     A
   |                  A A   |
   |          A  A  A       |
-18.338                       
   |                      |
-18.339            A    A     
   |                      |
   |   A                  |
-18.34                       
   |  A      A              |
-18.341      A               
   |                      |
       A               
-18.342 ----------------------------------------------------------
  1.78 1.782  1.784 1.786  1.788  1.79  1.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    |
            A  A   A     
-18.2 -------------------------------------BA--------------------
   1.7   1.72 1.74  1.76  1.78  1.8  1.82  1.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   10  12  14复制代码两端的数值会因偏离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
这个方法，即使在不高的水平下，结果也是比较精确的。

比较详细！

很详细，对初学者帮助很大。