Ex52 影响功函数计算的因素（二）

发表于 2018-09-06 更新于 2021-12-27 分类于 LVASPTHW 阅读次数： Disqus：

前面功函数的计算流程( VASP 计算细节以及结果分析)，我们已经掌握了。这一节介绍一下影响功函数的主要因素：真空层的厚度。首先个人经验：计算功函数的时候，真空层厚度一定要足够。

这里有两点：

多厚才算足够?具体要根据自己的体系测试下，一般来说 20 Å 左右足以。
是不是越厚越好? 越厚的话，计算出来的结果差别不大，但会增加计算量。前面我们讲过影响计算量的因素。

我们以 Cu(111) 表面的例子，简单介绍下真空层厚度对功函数的影响。

1 不同真空层厚度的测试(上一节内容)

先准备一个真空层为 10 Å 的计算文件夹： POSCAR 如上图；
将前面两节的 INCAR,KPOINTS, POTCAR, 以及提交任务的脚本复制过来；
批量生成测试文件，每个文件夹以真空层的厚度命名。命令如下：
1
for i in $(seq 12 2 36); do cp 10 $i; sed -i "5s/16/$((6+$i))/g"$i/POSCAR ; done
批量提交任务。

注意：

这样的批量操作本书刚刚开始的时候就已经讲到了，这里不要再问我。
如果不会用这个办法批量制作输入文件，想其他的办法批量操作，条条大路通罗马，不要死磕在我的这个命令上。

2 结果处理

前面的测试计算中，共有 14 个任务。在 Ex50 这一节中vtotav-v5.2f, vtotav.py 或者 p4vasp 这三种可视化的方法，我们查找真空能级的数值，都必须读取数据文件，比如 VLINE, LOCPOT_Z 或者 p4vasp 导出来的 *dat 文件。所以，我们要分别在这 14 个任务的文件夹中重复 Ex50 的操作14次。手动挨个弄可以完全实现，但本人不想这么做，因为这样弄除了让你的操作更熟练些外，并不会有其他的收获。那就写个脚本批量处理吧。

首先我们分析下上图：我们获取的真空能级在图中直线部分中间的一个点或者一个小区间的平均值。大约是真空层一半的高度。比如，我们的表层原子 z 方向坐标为 6，真空层厚度为 16，那么我们可以取 6 + (16/2)= 14 时对于的点来得到真空能级。也可以通过 [14-1, 14+1] 这个区间中所有的真空能级求平均数来获得。下面我们看一下脚本的主要内容：

注意部分：

使用这个脚本的时候，需要注意的内容。本人写了一堆废话在脚本里面，主要是提醒大家使用脚本的时候，需要注意的一些事项。

#!/usr/bin/env python

# Welcome to visit our website: www.bigbrosci.com to get more useful information.
 
'''
Written By Qiang 
This script has two functions:
 1) workfunction Visualization from LOCPOT_Z file
 2) Roughly estimate the Vaccum Energy

Please Contact lqcata@gmail.com if you have questions.

Note1: Read POSCAR with Cartesian coordinations! 
  This means 
  1) you have to convet the direct to cartesian firstly 
  2) otherwise you will get an error.

Note2: Read LOCPOT_Z to plot the figure and calculate the vaccum energy
  This means 
  1)1you have to run command:  work.py LOCPOT z firstly to get the LOCPOT_Z file before use this script

Note3: the idea to calculate the vaccum energy is following:
  1) get the middle z value in the plateau: 
     middle z value = ( coordination of the highest atom + length of the slab in z direction) / 2 
  2) selecte the area to calculate the vaccum energy: 
     from  (middle z  value - 1) to (middle z value + 1)  unit is in angstrom 
  3) do the average of all points y direction.

Note 4: this is only a rough estimation but useful.

Note 5: Check the figure firstly and then use the numbers calculated from this script. 
  1) If the midlle z value is far from the plateau in the figure, you have to calculate the energy by hand.

'''

脚本是谁写的，主要功能是什么；
脚本读取 POSCAR, POSCAR 要为 Cartesian 坐标；
脚本默认读取 LOCPOT_Z 文件，如果没有该文件，你需要做的一些事情(生成以及如何将 VLINE 文件转换为 LOCPOT_Z 文件)
该脚本是如何获取真空能级的；
一定要查看功函数可视化的结果，用来判断我们的这个方法是否试用。

写这么多废话，原因主要有 2 个：

提醒大家在使用别人脚本的时候，一些需要注意的事项你要清楚。如果不清楚，可以咨询下脚本的作者；
如果你写了一个脚本放到网上供大家使用，请花些时间写清楚脚本的运行原理，以及尽可能详细的注意事项。

脚本正文：

import matplotlib.pyplot as plt
import numpy as np 
x = []
y = []
dic = {}
with open("LOCPOT_Z", mode='r') as f:
    first_line = f.readline()
    name_x = first_line.split()[1]
    name_y = first_line.split()[2]
    for line in f:
        xy=line.rstrip().split()
        x.append(float(xy[0]))
        y.append(float(xy[1]))
        dic.update({xy[0]:xy[1]})
plt.plot(x,y) 
plt.xlabel(name_x) 
plt.ylabel(name_y) 
plt.savefig('workfunction' + '.pdf', dpi=400)
plt.show()

#%%%%%%%%%%%%%%%%%%%%%%#
# Get the Vaccum Energy
#%%%%%%%%%%%%%%%%%%%%%%#

# Get the total line numbers of POSCAR 
num_lines = sum(1 for line in open('POSCAR'))

# Read POSCAR 
pos = open('POSCAR', mode = 'r')
line = pos.readlines()

# Get  the  slab length in z direction 
vaccum = float(line[4].split()[2])

# Get all atoms' coordination in z direction and store them in the list
z_list = []
for i in range(9, num_lines): 
    z_list.append(float(line[i].split()[2]))

#  max(z_list) is highest atoms' coordination in z direction
# Get the vaccum lenth: 
l_vaccum = vaccum - max(z_list)
print 'The Vaccum in this calculation is:\t\t %s'  %(l_vaccum)

# Choose the middle z value in the workfuntion.pdf 
num_middle = (max(z_list) + vaccum) / 2 
#print num_middle

middle_y = []

for i in dic.keys():
    i = float(i)
# Select the date area within 1 angstrom from the middle point:
    if i > num_middle -1 and i <  num_middle + 1:
        middle_y.append(float(dic.get(str(i))))
# Get the average value in the selected area
print 'The Vaccum Energy in this calculation is:\t %s'  %(np.mean(middle_y))

pos.close()

plt.show() 之前为我们 Ex50 这一节可视化的脚本： 44 和 53 行中，额外将 LOCPOT_Z 中两列的数据保存到一个字典里面，便于后面根据 z 方向的坐标查找真空能级。

获取 POSCAR 文件总的行数

1	num_lines = sum(1 for line in open('POSCAR'))

读取 POSCAR 文件

1 2	pos = open('POSCAR', mode = 'r') line = pos.readlines()

获取 slab 在 z 方向的数值
1
vaccum = float(line[4].split()[2])

将所有原子的 z 坐标保存到 z_list 数列中

1
2
3

z_list = []
for i in range(9, num_lines): 
  z_list.append(float(line[i].split()[2]))

获得真空层的厚度
1
l_vaccum = vaccum - max(z_list)
计算距离表层原子一半真空层厚度的z方向的数值： num_middle
1
num_middle = (max(z_list) + vaccum) / 2

将 [num_middle-1, num_middle+1] 区间范围内所有的真空能级提取出来，并保存到middle_y 这个数列中

middle_y = []

for i in dic.keys():
    i = float(i)
# Select the date area within 1 angstrom from the middle point:
    if i > num_middle -1 and i <  num_middle + 1:
        middle_y.append(float(dic.get(str(i)))

求平均值，并打印出来

1	print 'The Vaccum Energy in this calculation is:\t %s' %(np.mean(middle_y))

运行脚本：

图中第一个命令中，输出比较复杂，并且会和前面示例一样，展示功函数可视化的结果。

但这不便于批量提取数据，可注销脚本中一下三行，直接获取真空能级
1
2
3
#plt.show()
#print 'The Vaccum in this calculation is:\t\t %s' %(l_vaccum)
#print 'The Vaccum Energy in this calculation is:\t %s' %(np.mean(middle_y))
使用脚本批量获取真空能级

先将脚本复制到bin文件夹中，赋予可执行权限。
1. 批量生成 LOCPOT_Z 文件：
2. 使用脚本批量读取
批量获得费米能级：

我们需要的为第 4 列的数据。(图中的冒号也算作一列)