计算材料学
0%

自动计算对称表面表面能脚本

发表于 更新于 分类于
本文字数： 2.5k 阅读时长 ≈ 2 分钟 
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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
# 20200422添加mp库修正原子能量
eles = {'La':-4.9266, 'Ce':-5.9343, 'Pr':-4.7837, 'Y':-6.4674,
        'Er':-4.5674, 'Li':-1.9080, 'Sm':-4.7163, 'Ni':-5.7783,
        'Cu':-4.0993, 'Ag':-2.8262, 'Sr':-1.6841, 'Eu':-10.2906,
        'Hg':-0.3016, 'Tm':-4.4782, 'Ho':-4.5781, 'Tb':-4.6275,
        'Lu':-4.5259, 'Dy':-4.5999, 'Ba':-1.9246, 'Nd':-4.7651, 
        'Pr':-4.7837, 'F':-1.882, 'Sc' : -6.3325,'H':-3.23975,
        'He': -0.0255} 

mu_mg = -1.594

import os

ele = str(input('请输入另一种元素符号(注意大小写格式):'))
index_a =eval(input('请输入化学式中Mg的原子个数：'))
index_b =eval(input('请输入化学式中另一种元素的原子个数：'))
formation_energy=eval(input('请输入化合物每原子的生成焓：'))
mu_formula_atom=eval(input('请输入化合物每原子的DFT能量：'))
print('----------计算结果输出如下----------')

'''
total_num:单胞中总共的原子个数
mu_formula：单胞DFT修正后能量
ratio：用来判断切好的表面中元素个数的关系，要减去多余的不符合化学计量比的原子
mu_F：在富Mg情况下第二种元素的原子化学势
low：假设的一个表面能初始值，1000.用于后面找到最小的表面
'''

total_num =index_a+index_b
mu_formula =mu_formula_atom*total_num
ratio = index_b/index_a
low = 1000.0
mu_F = eles[ele]+formation_energy*total_num/index_b

#获取当前目录
path = os.getcwd()
dirs = os.listdir(path)
#进入到表面目录进行循环
for dir in dirs:
    if dir[0]=='0' or dir[0]=='1':
        os.chdir(path+'/'+dir)
        temp_poscar = open("POSCAR", "r", encoding="utf-8").readlines()
        poscar = []
        for i in temp_poscar:
            poscar.append(i.strip())
        oszicar = open("OSZICAR", "r", encoding="utf-8").readlines()
        x=poscar[5][0:-1].split(" ")[0]
        y=poscar[5][0:-1].split(" ")[-1]
        #获取元素的系数，x是Mg
        if x[0] == 'M':
            num_x=eval(poscar[6].strip().split()[0])
            num_y=eval(poscar[6].strip().split()[1])
        else:
            num_x=eval(poscar[6].strip().split()[1])
            num_y=eval(poscar[6].strip().split()[0])
        #获取表面积
        temp_a = list(filter(None, poscar[2].split(" ")))
        temp_b = list(filter(None, poscar[3].split(" ")))
        area = eval(temp_a[0])*eval(temp_b[1])
        #获取收敛能量，但目前无法判断结构是否是收敛
        try:
            energy = eval(oszicar[-1].strip().split()[4])
        except:
            energy = eval(oszicar[-1].split(" ")[7])
        n = 0
        m = 0
        #判断系数关系
        if (num_y/num_x) == ratio:
            surf = (energy - num_x/index_a*mu_formula)/2/area
        elif (num_y/num_x) < ratio:
            while (num_y/num_x) < ratio:
                n = n+1
                num_x = num_x-1
            surf = (energy - num_x/index_a*mu_formula - n*mu_mg)/2/area
        else:
            while (num_y/num_x) > ratio:
                num_y = num_y-1
                m = m+1
            surf = (energy - num_x/index_a*mu_formula - m*mu_F)/2/area
        #输出表面能
        convert_surf = surf*16000
#         print(os.getcwd())
#         print(energy)
        print("{}面,{:.3f} mj/m2".format(os.getcwd().split('/')[-1],convert_surf))
        if low > convert_surf:
            low = convert_surf
        else:
            continue
    else:
        continue
# -0.018
# -2.4328
print('----------计算结果为----------')
print('表面能计算结束,最稳定表面为：'+str(round(low,3)))