1.3 结构优化3

发表于 更新于 分类于

加压

可通过施加静水压强、调整晶格常数的方法实现结构压缩

vi INCAR 增加标签 PSTRESS=50
（单位为 0.1 GPa）

其他不变

加电场

在vasp中根据自己所需要添加电场的场强更改INCAR中的如下参数：

1
2
3
LDIPOL=.TRUE.
IDIPOL=1/2/3/4
EFIELD= ？

对于IDIPOL参数的设置，需注意：对于IDIPOL=1-3，偶极矩将分别平行于第一、第二或第三晶格矢量的方向计算。总能量的修正计算为当前超单元中单极子/偶极子和四极子与放置在超单元中的同一偶极子之间的能量差，对应的晶格矢量接近无穷大。此标志应用于板计算。对于IDIPOL=4，将计算所有方向上的完整偶极矩，并将总能量的修正计算为当前超级单体中的单极/偶极/四极和放置在真空中的同一单极/偶极/四极之间的能量差，使用该标志计算孤立分子。

此外，EFIELD 的设置需要注意单位的换算。eV/Å（电子伏特/埃）计算完毕之后通过后处理软件分析所关注的性质即可。

状态方程拟合

我们需要对结构进行优化才可以获取稳定的晶格参数信息。有两个方法可以实现：

1 Birch-Murnaghan状态方程拟合,
2 VASP计算中通过调节ISIF参数直接优化Bulk。
直接优化前面已经介绍，下面我们先讨论一下第一个方法：BM方程拟合。

参考：大师兄科研网：https://www.bigbrosci.com/2018/02/02/ex33/

科学网耿华运的博文：https://blog.sciencenet.cn/blog-3469498-1280203.html

科学网赵志强的博文，关于拟合：https://wap.sciencenet.cn/blog-3388193-1152915.html

vaspkit：http://vaspkit.cn/index.php/48.html

大致流程（以大师兄提供脚本为例）

  1. 获得一定缩放系数的结构：for i in 0.95 0.96 0.97 0.98 0.99 1.01 1.02 1.03 1.04 1.05; do cp 1.00 $i ; sed -i "2s/1.0/$i/g" $i/POSCAR ; done

  2. 分别进行结构优化:for i in *; do cd $i; 提交任务命令; cd $OLDPWD; done

  3. 提取能量-体积以及能量-压力、压力-体积的数据:for i in *; do echo -e $i "\t" $(grep ' without' $i/OUTCAR | tail -n 1| awk '{print $7}'); done > data

  4. 拟合能量体积曲线，获得平衡晶胞参数

  • 这个流程可通过批量计算脚本实现，参考上面链接。另外还可以获得能量随压力的变化来模拟相变。

分数占据

生成disorder或者掺杂构型的软件汇总：

https://mp.weixin.qq.com/s/vYar9UcYWAZgazBTgqUGEA
https://github.com/jichunlian/disorder
supercell： https://orex.github.io/supercell/
ATAT： https://www.brown.edu/Departments/Engineering/Labs/avdw/atat/
vaspkit：https://vaspkit.com/

disorder

安装

1
2
3
4
5
6
7
8
$ unzip disorder-main
$ cd disorder-main
$ make disorder (Only disorder is compilied)
$ make supercell (Only supercell is compilied)
$ make or make all (Both disorder and supercell are compilied)
$ make clean (rm *.o *.mod disorder supercell)
export PATH=$PATH:disorder_installation_path/disorder-main/bin
(Add this to the ~/.bashrc file as you wish)

Note: The default compiler is ifort, and you can modify the Makefile
in the src directory to use another compiler (e.g. gfortran).

需要准备两个文件SPOSCAR和INDSOD,以examples的1和2为例：

例子1：SnxPb1-xTe

  1. 扩建超胞，保证整除
    1
    2
    3
    4
    5
    6
    7
    8
    9
    Te Pb
       1.0
        13.1323995589999996    0.0000000000000000    0.0000000000000000
         0.0000000000000000   13.1323995589999996    0.0000000000000000
         0.0000000000000000    0.0000000000000000   13.1323995589999996
      Te   Pb
      32   32
    Direct
    * * *
  2. 准备INDSOD文件
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    &input
      nsub = 2
      subs = 16,16
      symb = 'Sn','Pb' ! The quotes is unnecessary for the ifort compiler
      site = 2
      prec = 1D-5
    ! lpro = .true.
    ! lpos = .true.
    ! leqa = .true.
    ! lspg = .true.
      lcfg = .false.
    /
    官方解释：
    1
    2
    3
    4
    5
    6
    7
    nsub-tag:nsub = 2 ~ 5; Default: nsub = 2.该位置取代原子种类。本例中是2种。
    subs-tag：subs = integer, integer, …; Default: None。对 应的原子数目。总和加起来应该等于该位点在SPOSCAR中的数目。
    symb-tag：symb = character, character, …; Default: None。元素化学符号。空位用Kw,即汉语kong wei。本例是Sn Pb单引号。
    site-tag：site = integer; Default: 1。所取代的位置在SPOSCAR中的位置。本例是第2.
    prec-tag：prec = real; Default: 1D-5。当软件产生空间群与其他软件识别不一致时调整该参数，一般不用修改。
    lpro-tag：lpro = logical; Default: .false.是否消除重复参数。对于大体系可以打开。
    lpos、leqa、lspg、lcfg：分别控制输出文件POSCAR-X、EQAMAT、SPGMAT、CFGMAT。
  3. 运行disorder。输入disorder即可。
  4. 输出文件：

ATAT软件

简介

  • 什么是ATAT？

其全称为 Alloy Theoretic Automated Toolkit (ATAT)
引用官网的一句话，ATAT is a generic name that refers to a collection of alloy theory tools developped by Axel van de Walle, in collaboration with various research groups and with various sources of financial support.

  • 如何安装ATAT：

http://brown.edu/Departments/Engineering/Labs/avdw//atat/atat3_36.tar.gz下载ATAT源码，解压后，make all开始编译；之后make install进行程序和脚本定向安装；之后在home中打开./bashrc添加ATAT安装路径于环境变量当中（添加 export PATH=$PATH:/xxxx），并source ./bashrc即可。

进入安装包的glue目录，里面有与各种dft软件的算例的接口设置。例如进入vasp可以看到它的接口设置脚本ezvasp，由于make install已经将它拷贝到了安装目录中，所以在bin中打开ezvasp，并设置相关可执行程序，包括赝势文件的名称和路径即可。

配置~/.ezvasp.rc

1
2
3
#!/bin/csh
set VASPCMD="/path/vasp_std# 设置vasp可执行程序
set POTLDA/PBE="/path/to/pot_lad" # 设置为赝势实际的路径
  • ATAT能做啥？

ATAT 是Walle及其合作者开发的一个利用第一性原理结合MonteCarlo模拟计算相图及热力学性能的软件包。该软件包首先采用团簇展开(cluster expansion)，用第一原理方法计算一系列构型的形成能，然后用这些能量作为Monte Carlo模拟的哈密顿量，计算材料的热力学性能和高熵合金的mcsqs建模。

MAPs

软件运行

一）运行MAPS：

  1. 准备lat.in文件

2）maps -d & #maps命令运行，等待后面的指令

3）touch ready

#等待十几秒钟，出现“0”文件夹，是合金中比例为0的纯相结构。

4）此时需要利用第一性原理软件计算该结构的能量

#（参见下部分和VASP软件的交互,可无需手动输入），现在手动输入能量，理解一下计算过程，假设该结构的能量已经知道，在0这个文件夹下，echo 1.1 > energy, 并删除wait。继续touch ready，就会生成1这个文件夹，里面是比例为1的另一个纯相结构，echo 2.2 > energy, 删除wait。touch ready，生成2这个文件夹，一直这样循环下去，直至达到相应的精度，最终停止输入touch stop。

二）运行MAPS与VASP交互：

1）需要lat.in，vasp.wrap文件 （ #vasp的INCAR参数）

2）其余与上面运行MAPS步骤相同，直到第4步骤计算能量，进入该文件夹，运行runstruct_vasp。

大概率会报错，原因有几个：

a. 常见为~/.ezvasp.rc文件配置问题，赝势位置，vasp运行路径出错。

可以看运行命令后，是否有POTCAR文件产生，如果没有，基本就是赝势位置问题，如果POTCAR文件正常，就是vasp没有运行。

b. vasp不运行，在集群上面运行的软件，ATAT软件给了集群上面的运行办法, pollmach命令，将runstruct_vasp命令写在脚本里是可以的。

最终正确版本运行步骤：

1）准备lat.in，vasp.wrap文件，脚本文件。

lat.in文件

vasp.wrap文件

sub.sh文件

2）提交脚本文件即可

说明：

  • 缺少.machines.rc文件的提示可忽略，该文件用于跨节点的并行
  • maps不断产生结构和构建团簇展开模型，-d表示使用默认参数，直接输入maps可查看可用参数
  • runstruct_vasp调用vasp计算所产生的结构的能量，每个目录包含一个结构及其计算结果，每个目录下的energy文件即从vasp的计算中读取的能量
  • maps.log中记载了现有的团簇展开模型的状态，包括对基态构型的预测和LOOCV等

什么时候停止计算：

  • 团簇展开模型能够成功预测基态 (查看maps.log)
  • LOOCV足够小 (<0.025 eV suggested by ATAT manual)
  • 有效团簇相互作用随着团簇距离的增大而衰减

以上的三个判据并不是绝对的，判断一个团簇展开模型是合理的还可以有其它方式。

停止计算：

1
2
$ touch stop
$ touch stoppoll

重启计算: 与正常开始一个计算相同，如果之前的计算不是正常结束则需要删除一些文件

1
2
3
rm maps_is_running pollmach_is_running ready
maps -d &
pollmach runstruct_vasp

常用命令

  • corrdump -2= -3= -4= 产生团簇

  • getclus 查看团簇

    • 第一列为团簇的阶数，第二列为团簇直径，第三列为团簇的多重度
    • clusters.out中有团簇格点位置的详细信息

  • clusterexpand -e -cv energy 得到团簇展开系数

    • 此时产生的ECI文件为energy.eci，而maps产生的ECI为eci.out，后者是形成能的团簇展开结果，前者是直接对vasp能量(即energy文件)的团簇展开
    • 如果需要对其它性质做团簇展开，例如带隙，则可以在每个结构的目录下创建一个bandgap文件,其内容即为带隙的大小，然后运行clusterexpand -e -cv -pa bandgap以得到带隙团簇展开的系数bandgap.eci，其中-pa表示所展开的量已经是一个平均的量

检查计算结果

  • VASP是否正常结束
    • ATAT会自动检查一些错误(checkerr_vasp)
    • 但如果结构优化的离子步数不够大，不会有提示，因此建议在vasp.wrap中将NSW设置大一些
  • 晶胞的应变是否很大
    • checkrelax (<0.1 suggested by ATAT）

结果可视化

  • mapsrep 会调用gnuplot作图并产生相应的脚本mapsrep.gnu

输出文件

  1. maps.log：开始可能的warnings:
    ’Not enough known energies to fit CE’
    ’True ground states not = fitted ground states’
    ’New ground states predicted, see predstr.out’
    随着结构增多:
    ’Among structures of known energy, true and predicted ground states agree.’
    ’No other ground states of xx atoms/unit cell or less exist.’
    最后一行有crossvalidation score (in eV/atom)，该值要小于0.025。精确小于0.001.

  2. fit.out:拟合的结构，一个结构一行，共6列，分别代表concentration、 energy 、fitted_energy (energy-fitted_energy)、 weight 、index

  3. predstr.out：预测的能量，一个结构一行，共5列，分别代表concentration、 energy 、predicted_energy、index、status.其中status包括status is either
    b for busy (being calculated),
    e for error,
    u for unknown (not yet calculated) or
    g if that structure is predicted to be a ground state (g can be combined with the above).

  4. gs.out：基态能量，一个结构一行，四列，分别是concentration、 energy、 fitted_energy、 index

  5. gs_str.out：基态的结构，格式和str.out一样。每一个结构后面会有end和一个空行隔开。

  6. eci.out：eci值对应的cluster在clusters.out文件。

  7. clusters.out：对于每一个团簇，第一行是多重度 multiplicity，第二行是团簇维度cluster diameter，第三行是团簇的格点数，下一行是格点对应的坐标。最后空行隔开每一个团簇。

  8. str.out：在每个数字文件夹下面，和晶格输入文件相似，不同是坐标写成3*3矩阵，一个位置只写一个原子。

  9. ref_energy.out：用于计算形成能的参考能量。

后处理

一系列用于分析ATAT计算结果的python脚本，GitHub地址为：https://github.com/xu-xi/clusterexpansion可以通过 git clone 下载到本地服务器。用“-h”命令可以查看每个脚本的帮助文档。最常用的几个脚本是：

1）cv_scan.py，测试CE模型对团簇截断半径的收敛性。

2）lc.py，测试CE模型对训练集大小的收敛性。通常~80个构型就可以得到比较精确的团簇展开模型。

3）ce_plot.py，这个脚本绘制的是拟合能量-计算能量分布图，默认-p=energy（调用的是每个数字目录下的energy文件）。energy文件是ATAT完成一个构型的DFT (relax+static) 计算后自动生成的能量文件。在实际应用中，拟合能量-计算能量分布图需要的往往不是energy而是formation energy，所以就需要编写脚本手动在每一个数字目录下生成一个形成能文件formation_energy，该值需要将目标目录下energy的值与0文件夹的energy和1文件夹的energy作线性变换得到，再写入到每一个数字目录中，读者可自行尝试。下面的脚本仅供参考，

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#!/usr/bin/env python
import os,subprocess,numpy,argparse,sys,subprocess,math
energy_0 = float(subprocess.check_output('cat ./0/energy', shell = True))
energy_1 = float(subprocess.check_output('cat ./1/energy', shell = True))
for item in os.listdir(os.getcwd()):
    fullpath = os.path.join(os.getcwd(),item)
    if os.path.isdir(fullpath):
        os.chdir(fullpath)
        if not os.path.isfile('error'):
            sc_atom_number = len(file('str.out').readlines())-6
            # print sc_atom_number
            energy = float(subprocess.check_output('cat energy', shell = True))
            energy /= sc_atom_number
            #energy *= 5 normalize to unit cell,5 atoms in CsPbI3
            sc_x = float(subprocess.check_output('corrdump -pc -l=../lat.in',shell=True).strip().split()[0])
            formation_energy = energy - (1 - sc_x) * energy_0 - sc_x * energy_1
            # print formation_energy
            with open('formation energy','w') as f:
                f.write('%.7f\n' % formation_energy)
        os.chdir('../')

注：这里如果报错：NameError: name ‘file’ is not defined. Did you mean: ‘filter’?将file改成open

得到formation_energy文件后，执行命令ce_plot.py -p=formation_energy -s=”formation energy”，还可以用Origin手动绘制出形成能-掺杂浓度分布图；在得到了比较精确的CE模型后，可以进行蒙特卡罗模拟，得到团簇关联函数、系综平均能量和方差随温度变化示意图。

热力学计算

必要的输入文件：

  • lat.in
  • clusters.out
  • eci.out
  • gs_str.out
emc2使用示例

蒙特卡洛模拟

1
emc2 -T0=300 -T1=3000 -dT=100 -cm -x=0 -keV -gs=-1 -er=30 -aq=0 -dx=1e-5 -sigdig=12
  • -cm表示正则系综
  • -x其实是点团簇函数，若c表示格点的占据浓度，则x=2c-1
  • -sigdig=12是为了表示有足够的位数输出能量的方差以判断相变
  • -aq和-dx是自动确定MC平衡步数和取样步数的算法，其原理可参考文献Ref.1，也可以用-eq和-n来直接指定步数，需要做收敛性测试
  • 主要输出文件为mc.out，输入emc2 -h可查看每一列所表示的内容
  • 输出文件mcsnapshot.out为最后一个模拟温度下的MC快照。如果想得到某个温度下的snapshot，可以只设置T0而不设置T1和dT。 
    1
    2
    3
    $ cp mcsnapshot.out str.out
    $ runstruct_vasp -nr 
    $ vesta POSCAR

相图计算

基本原理可参考文献Ref.1，Martin Baker的笔记提供了更多的例子和细节

例子：

1
phb -dT=10 -ltep=1e-3 -er=30 -gs1=0 -gs2=1 -o=phb.out -keV -dx=1e-5

构建SQS/SQoS

SQS/SQoS是对无序体系的代表性结构，其中SQS假定原子完全随机的占据，而SQoS是SQS的扩展，可以考虑一定的短程有序性，其基本原理可参考相关文献Ref.2

必要文件：

  • rndstr.in
  • clusters.out

rndstr.in的一个例子：

1
2
3
4
5
6
7
8
9
10
11
4.1000 0.000000 0.000000
0.000000 4.1000 0.000000
0.000000 0.000000 4.1000
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
0.000000 0.500000 0.500000 O=0.6667,N=0.3333
0.500000 0.000000 0.500000 O=0.6667,N=0.3333
0.500000 0.500000 0.000000 O=0.6667,N=0.3333
0.000000 0.000000 0.000000 Ba
0.500000 0.500000 0.500000 Ta

产生团簇文件clusters.out：

1
corrdump -l=rndstr.in -clus -noe -nop -2= -3= -4=

说明：

  • -noe和-nop表示忽略空团簇和点团簇
  • 用getclus查看所产生的团簇。如果团簇数目过少，则后续产生SQS有可能很容易完全匹配团簇函数，导致SQS实际上不够随机；而如果团簇过多，则有可能导致近邻的团簇匹配得不够好，而通常近邻的团簇是更为重要的，可以通过-wd加上权重以避免这种情况

运行mcsqs -n= &
说明：

  • -n设定SQS晶胞的原子数
  • 可通过-tcf设定目标团簇函数（可从mc.out中提取），即SQoS方法
  • mcsqs不会自动停下来，可用touch stopsqs停止
  • 输出文件mcsqs.log中有对团簇函数匹配的优化记录
  • 输出文件bestsqs.out为目前找到的最佳的SQS，可将其复制为str.out，然后运行runstruct_vasp -nr生成vasp的输入文件用于后续计算

MCSQS

一、原理

无序态的关联函数已知，因此只需要备选态的关联函数和无序态关联函数差别最小就可以，即让下方公式中的目标函数最小。

二、运行步骤

step1. 准备rndstr.in文件

1 1 1 90 90 90

1 0 0

0 1 0

0 0 1

.0 .0 .0 Cr

.0 .5 .5 Ni=0.333333,Fe=.333333, Co=.333334

.5 .0 .5 Ni=0.333333,Fe=.333333, Co=.333334

.5 .5 .0 Ni=0.333333,Fe=.333333, Co=.333334

#第一行也可以改成3*3矩阵，类似于POSCAR中晶格常数格式。这里一般不用实际的晶格常数，用1替代。

用1对后面的参数设置更容易。

step2. 运行命令
corrdump -l=rndstr.in -ro -noe -nop -clus -2=1.1
该命令通过读取rndstr.in文件中的结构，进而计算结构的对称性和团簇等信息，输出clusters.out文件。

#需要注意的是‘-2’这个参数，是用于计算关联函数的截断距离。如对于FCC构型而言，晶格常数为1，第一最近距离为2^0.5/2=0.71，第二距离为1，第三距离为1.5^0.5=1.2。所以这里设置阶段距离在1和1.2之间。一般考虑‘-2’参数即可。（注：更高的-2，-3参数会产生更无序的结构，但是需要更长的时间，一般晶格常数为1，这里选取1.1就可以）

step3. 运行命令 mcsqs -n XX
该命令主要输出sqscell.out文件。XX表示体系的原子个数，该数值一定是上述结构原子的倍数，否则会出错。

#比如上述构型，自身CrNi有4个原子，其中Ni元素与Fe,Co构成合金，因此XX至少为12个原子。该步骤运行一下，就可以进行停止。输出文件需要关注sqscell.out文件，该文件给出了XX个原子扩胞后形成的不同晶格构型。可通过挑选，删除一些不太可能的结构，比如x,y,z三个方向长度差别过大，形成的角度较小等对VASP或后续计算不友好的结构。之后用新的sqscell.out文件运行第四步骤。

step4. 运行命令mcsqs -rc
输出bestsqs.out文件和bestcorr.out文件。其中bestsqs.out就是无序结构。

#通常而言没有最完美的无序结构，这个命令会一直运行下去，可以肯定的是，再一次生成的结构比上一次生成的结构要优异，所以什么时候停止取决于自身研究的精度。

层间结合能、结合力

贴一篇文章：https://news.sciencenet.cn/htmlpaper/2022/4/202242914585844872440.shtm

Graphite interlayer distance

vasp上的教程：

准备graphite结构，

INCAR：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
IVDW = 20           
LVDW_EWALD =.TRUE. 
NSW = 1 
IBRION = 2
ISIF = 4
PREC = Accurate
EDIFFG = 1e-5
LWAVE = .FALSE.
LCHARG = .FALSE.
ISMEAR = -5
SIGMA = 0.01
EDIFF = 1e-6
ALGO = Fast
NPAR = 2

提交脚本：
To run this example, execute the run.sh bash-script:

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
#
# To run VASP this script calls $vasp_std
# (or posibly $vasp_gam and/or $vasp_ncl).
# These variables can be defined by sourcing vaspcmd
. vaspcmd 2> /dev/null

#
# When vaspcmd is not available and $vasp_std,
# $vasp_gam, and/or $vasp_ncl are not set as environment
# variables, you can specify them here
[ -z "`echo $vasp_std`" ] && vasp_std="mpirun -np 8 /path-to-your-vasp/vasp_std"
[ -z "`echo $vasp_gam`" ] && vasp_gam="mpirun -np 8 /path-to-your-vasp/vasp_gam"
[ -z "`echo $vasp_ncl`" ] && vasp_ncl="mpirun -np 8 /path-to-your-vasp/vasp_ncl"

#
# The real work starts here
#

rm results.dat

for d in 6.5 6.6 6.65 6.7 6.75 6.8 6.9 7.0
do
  cat>POSCAR<<!
graphite
1.0
1.22800000 -2.12695839  0.00000000
1.22800000  2.12695839  0.00000000
0.00000000  0.00000000  $d
4
direct
   0.00000000  0.00000000  0.25000000
   0.00000000  0.00000000  0.75000000
   0.33333333  0.66666667  0.25000000
   0.66666667  0.33333333  0.75000000
!
  $vasp_std
  cp OUTCAR OUTCAR.$d
  energy=$(grep "free  ene" OUTCAR.$d|awk '{print $5}')
  echo $d $energy >> results.dat 
done

The optimal length of the lattice vector c normal to the stacking direction is determined in a series of single point calculations with varied value of c (all other degrees of freedom are fixed at their experimental values).

The computed c vs. energy dependence is written in the file results.dat and can be visualized e.g. using xmgrace. The optimal value can be obtained using the attached utility (python with numpy or Numeric is needed):

./utilities/fit.py results.dat
This should yield:

1
2
3
4
5
6
7
200 iterations performed
Ch-square: 4.30305519481e-09
---------

       E0(eV):         -37.433456779
       d0(A):  6.65603352689
The computed value of 6.66 Å agrees well with experiment (6.71 Å).

Graphite MBD binding energy

1
2
3
4
5
6
7
8
9
10
11
12
13
14
IVDW = 202           
LVDWEXPANSION =.TRUE. 
NSW = 1 
IBRION = 2
ISIF = 4
PREC = Accurate
EDIFFG = 1e-5
LWAVE = .FALSE.
LCHARG = .FALSE.
ISMEAR = -5
SIGMA = 0.01
EDIFF = 1e-6
ALGO = Fast
NPAR = 2

KPOINTS

Graphite:

1
2
3
4
5
Monkhorst Pack
0
gamma
16 16 8
0 0 0

Graphene:

1
2
3
4
5
Monkhorst Pack
0
gamma
16 16 1
0 0 0

Running this example
To run this example, execute the run.sh bash-script:

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
#
# To run VASP this script calls $vasp_std
# (or posibly $vasp_gam and/or $vasp_ncl).
# These variables can be defined by sourcing vaspcmd
. vaspcmd 2> /dev/null

#
# When vaspcmd is not available and $vasp_std,
# $vasp_gam, and/or $vasp_ncl are not set as environment
# variables, you can specify them here
[ -z "`echo $vasp_std`" ] && vasp_std="mpirun -np 8 /path-to-your-vasp/vasp_std"
[ -z "`echo $vasp_gam`" ] && vasp_gam="mpirun -np 8 /path-to-your-vasp/vasp_gam"
[ -z "`echo $vasp_ncl`" ] && vasp_ncl="mpirun -np 8 /path-to-your-vasp/vasp_ncl"

#
# The real work starts here
#

# Here the work starts
rm results.dat

drct=$(pwd)

for i in graphene graphite
do
  cd $drct/$i
  $vasp_std
done

cd $drct

# obtain total energy for graphite 
en2=$(grep "free  ene" graphite/OUTCAR |tail -1|awk '{print $5}') 

# obtain total energy for graphene
en1=$(grep "free  ene" graphene/OUTCAR |tail -1|awk '{print $5}')

# compute interlayer binding energy (eV/atom)
deltaE=$(echo print $en2/4 - $en1/2 |python)

echo "Binding energy (eV/atom): " $deltaE >results.dat

Note that the calculation is performed in two steps (two separate single-point calculations) in which the energy for bulk graphite and for graphene are obtained. The binding energy is computed automatically and it is written in the file results.dat. (N.B.: for the latter python needs to be available.)

The computed value of 0.050 eV/A is now fairly close to the RPA reference of 0.048 eV/atom [1].

Graphite TS binding energy

1
2
3
4
5
6
7
8
9
10
11
12
13
14
IVDW = 20            
LVDW_EWALD =.TRUE. 
NSW = 1 
IBRION = 2
ISIF = 4
PREC = Accurate
EDIFFG = 1e-5
LWAVE = .FALSE.
LCHARG = .FALSE.
ISMEAR = -5
SIGMA = 0.01
EDIFF = 1e-6
ALGO = Fast
NPAR = 2

To run this example, execute the run.sh bash-script:

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
#
# To run VASP this script calls $vasp_std
# (or posibly $vasp_gam and/or $vasp_ncl).
# These variables can be defined by sourcing vaspcmd
. vaspcmd 2> /dev/null

#
# When vaspcmd is not available and $vasp_std,
# $vasp_gam, and/or $vasp_ncl are not set as environment
# variables, you can specify them here
[ -z "`echo $vasp_std`" ] && vasp_std="mpirun -np 8 /path-to-your-vasp/vasp_std"
[ -z "`echo $vasp_gam`" ] && vasp_gam="mpirun -np 8 /path-to-your-vasp/vasp_gam"
[ -z "`echo $vasp_ncl`" ] && vasp_ncl="mpirun -np 8 /path-to-your-vasp/vasp_ncl"

#
# The real work starts here
#

rm results.dat

drct=$(pwd)

for i in graphene graphite
do
  cd $drct/$i
  $vasp_std
done

cd $drct

# obtain total energy for graphite 
en2=$(grep "free  ene" graphite/OUTCAR |tail -1|awk '{print $5}') 

# obtain total energy for graphene
en1=$(grep "free  ene" graphene/OUTCAR |tail -1|awk '{print $5}')

# compute interlayer binding energy (eV/atom)
deltaE=$(echo print $en2/4 - $en1/2 |python)

echo "Binding energy (eV/atom): " $deltaE >results.dat

Note that the calculation is performed in two steps (two separate single-point calculations) in which the energy for bulk graphite and for graphene are obtained. The binding energy is computed automatically and it is written in the file results.dat. (N.B.: for the latter python needs to be available.)

Even though the TS method predicts a reasonable geometry (see the Graphite interlayer distance example) it overestimates the energetics strongly: the computed binding energy of -0.083 eV/atom is too large compared to the RPA reference of 0.048 eV/atom. This overestimation is - at least in part - due to neglecting the many-body interactions (see example Graphite MBD binding energy).