以p(1x1)Cu(111)为例分析表面弛豫的结果

本文其实没啥新的内容，只是对上两节的一个总结分析

在表面弛豫前（slab模型优化计算前）

1
2
3
4
5
6
7
8
9
10
11
12
13
SYSTEM = p(1x1)Cu(111)
1.0
        2.5701000690         0.0000000000         0.0000000000
       -1.2850500345         2.2257719501         0.0000000000
        0.0000000000         0.0000000000        21.2954006195
   Cu
    4
Selective Dynamics
Direct
     0.000000000         0.000000000         0.000000000 F F F
     0.666670000         0.333330000         0.098540000 F F F
     0.333330000         0.666670000         0.197080000 T T T
     0.000000000         0.000000000         0.295620000 T T T

优化计算后

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
SYSTEM = p(1x1)Cu(111)
1.0
        2.5701000690         0.0000000000         0.0000000000
       -1.2850500345         2.2257719501         0.0000000000
        0.0000000000         0.0000000000        21.2954006195
   Cu
    4
Selective Dynamics
Direct
     0.000000000         0.000000000         0.000000000 F F F
     0.666670000         0.333330000         0.098540000 F F F
     0.333330000         0.666670000         0.197080000 T T T
     0.000000000         0.000000000         0.295620000 T T T
[ctan@baifq-hpc141 primitive-slab-opt-1]$ cat CONTCAR
CONTCAR\(1\1\1)
   1.00000000000000
     2.5701000690000000    0.0000000000000000    0.0000000000000000
    -1.2850500345000000    2.2257719501000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.2954006195000005
  Cu/0e71e558f37
               4
Selective dynamics
Direct
  0.0000000000000000  0.0000000000000000  0.0000000000000000   F   F   F
  0.6666700000000034  0.3333299999999966  0.0985399999999998   F   F   F
  0.3333311322281505  0.6666688677718494  0.1972481653092927   T   T   T
 -0.0000004145140272  0.0000004145140272  0.2951613337135311   T   T   T

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00

工作站离线状态下numpy on python2地安装

这两个文件中的坐标都是分数坐标，不太利于分析位置变化，需要将分数坐标转化为笛卡尔坐标，使用大师兄的脚本来进行转换，

用的是python2

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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#Convert direc coordiation to cartesian Writen By Qiang and Xu Nan
"""
The right input way to run this script is : 
      
python dire2cart.py CONTCAR 3 

CONTCAR will be converted from Direct to Cartesian and 3 layers will be fixed.  

If you do not want to fix layers and keep the same as before, run command:
    
python dire2cat.py CONTCAR    

"""

import sys
import os
import numpy as np 

print 'Please read the head part of this script and get more information!'
print """
###################################
#                                 #
#for VASP 5.2 or higher versions  #
#                                 #
###################################
"""

if len(sys.argv) <= 1:
    print('\n' + ' Warning! ' * 3 + '\n')
    print('You did not select the inputfile to be converted. \n By defalut, we are going to convert your CONTCAR.\n')
    
    if not os.path.isfile("POSCAR") and not os.path.isfile("CONTCAR"):
        print("Error:" * 3 + "\n Can not find neither POSCAR nor CONTCAR!!\n")
        exit()
    else:
        if os.path.isfile("CONTCAR"):
            script = sys.argv[0]
            file_to_be_converted = "CONTCAR"
            print "\n Convertion starts........\n"
    
        else:
            script = sys.argv[0]
            file_to_be_converted = "POSCAR"
            print "\n There is no CONTCAR in your directory. \n \n POSCAR Convertion starts........\n"
            
if len(sys.argv) == 2:
    print("\n%s Convertion starts......" %sys.argv[1]) 
    script, file_to_be_converted = sys.argv[:2]
else:
    print("\n%s Convertion starts......\n" %sys.argv[1]) 
    script, file_to_be_converted, fixedlayer = sys.argv[:3]
    fixedlayer = int(fixedlayer)

def get_infor():
    f = open(file_to_be_converted, 'r')
    lines = f.readlines()
    f.close()
    num_atoms = sum([int(x) for x in lines[6].split()])
    if lines[7][0]  == 'S' or lines[7][0]  == 's':  # # With Selected T T T, coordination starts from line 9
        start_num = 9
        if lines[8][0]  == 'D' or lines[8][0]  == 'd':
            is_direct = True 
        else:
            is_direct = False
    else: 
        start_num = 8
        print '----------------------------------------------------'
        print 'Pay Attetion! There is no TTT in  %s   ' %(file_to_be_converted)
        print '----------------------------------------------------'
        if  lines[7][0]  == 'D' or lines[7][0]  == 'd':
            is_direct = True 
        else:
            is_direct = False
    a = []
    b = []
    c = []
    if is_direct:         
        for i in np.arange(2,5):
            line = [float(i) for i in lines[i].strip().split()]
            a.append(line[0])
            b.append(line[1])
            c.append(line[2])
        vector = np.array([a,b,c])
    else: 
        vector = np.array([[1, 0 , 0], [0, 1, 0], [0, 0, 1]])
    return vector, lines, start_num, num_atoms, is_direct

def determinelayers(z_cartesian):
    threshold = 0.5  
    seq = sorted(z_cartesian)
    min = seq[0]
    layerscount = {}
    sets = [min]
    for j in range(len(seq)):
        if abs(seq[j]-min) >= threshold:
            min = seq[j]
            sets.append(min)

    for i in range(1,len(sets)+1):
        layerscount[i] = []            
        for k in range(len(z_cartesian)):   
            if abs(z_cartesian[k]-sets[i-1]) <= threshold:
                layerscount[i].append(k)

    return layerscount

def convert():
    x_cartesian = []
    y_cartesian = []
    z_cartesian = []
    tf = []
    for i in range(start_num, num_atoms + start_num):
        line_data =  [float(ele) for ele in lines[i].split()[0:3]]
        line_data = np.array([line_data])
        x, y, z = [sum(k) for k in line_data * vector ]     
        x_cartesian.append(x)
        y_cartesian.append(y)
        z_cartesian.append(z)
         
        if start_num == 9 : # if  T T T exist, the start_num will be 9        
            tf.append((lines[i].split()[3:]))
        else:
            tf.append(' ')   # if there is no T T T, use space instead. 

    layerscount =determinelayers(z_cartesian)
    
    file_out = open(file_to_be_converted+'_C', 'w')
    for i in range(0,7):
        file_out.write(lines[i].rstrip() + '\n')  # first 7 lines are kept the same 
    
    print '\n Find %s layers!' * 3 %(len(layerscount), len(layerscount), len(layerscount))
    
    if len(sys.argv) >= 3:  # This means that the nuber for fixing layers is given by the user.
        file_out.write('Selective\n')  
        file_out.write('Cartesian' + '\n') 
        for i in range(1,len(layerscount)+1):
            if i <= fixedlayer: 
                for j in layerscount[i]:
                    tf[j] = ['F','F','F']
            else:
                for k in layerscount[i]:
                    tf[k] = ['T','T','T']                    
    else:
        if start_num == 9:  # 9 means there are T T T or F F F in the file
            file_out.write('Selective\n')
            file_out.write('Cartesian' + '\n')        
        else:
            file_out.write('Cartesian' + '\n') 
    
    for i in range(0,len(x_cartesian)):
        file_out.write("\t%+-3.10f   %+-3.10f   %+-3.10f  %s\n" %(x_cartesian[i], y_cartesian[i], z_cartesian[i], ' '.join(tf[i])))
        
    file_out.close()
    
vector, lines, start_num, num_atoms, is_direct = get_infor()

if is_direct : 
    print "\n%s has Direct Coordinates, Contersion starts.... "  %(file_to_be_converted)
    convert()
else:
    print "\n%s has Cartesian Coordinates Already! We are going to fix layers only."  %(file_to_be_converted)
    convert()
    
print '-----------------------------------------------------\n'
print '\n %s with Cartesian Coordiates is named as  %s_C\n' %(file_to_be_converted, file_to_be_converted)
print '-----------------------------------------------------\n'

脚本会调用到numpy库，然而集群上并没有连接互联网，所以没有办法通过pip2 install numpy来安装，只能手动上传.whl文件或者源码编译等。

进入python2环境，输入import pip;print(pip.pep425tags.get_supported())来查看支持的.whl文件版本

1
2
3
4
5
6
[ctan@baifq-hpc141 primitive-slab-opt-1]$ python2
Python 2.7.18 (default, Oct 23 2023, 20:01:03)
[GCC 8.5.0 20210514 (Red Hat 8.5.0-18)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import pip;print(pip.pep425tags.get_supported())
[('cp27', 'cp27mu', 'manylinux1_x86_64'), ('cp27', 'cp27mu', 'linux_x86_64'), ('cp27', 'none', 'manylinux1_x86_64'), ('cp27', 'none', 'linux_x86_64'), ('py2', 'none', 'manylinux1_x86_64'), ('py2', 'none', 'linux_x86_64'), ('cp27', 'none', 'any'), ('cp2', 'none', 'any'), ('py27', 'none', 'any'), ('py2', 'none', 'any'), ('py26', 'none', 'any'), ('py25', 'none', 'any'), ('py24', 'none', 'any'), ('py23', 'none', 'any'), ('py22', 'none', 'any'), ('py21', 'none', 'any'), ('py20', 'none', 'any')]

这里下载的是numpy-1.16.5-cp27-cp27mu-manylinux1_x86_64.whl（不可下载过高版本的numpy，它支持python2的最后一个版本似乎是1.16.6）

然后安装，由于我不在sudoer中，没有root权限，所以需要手动指定安装目录

1
pip2 install --prefix=~/python2/numpy_installation numpy-1.16.5-cp27-cp27mu-manylinux1_x86_64.whl

然后添加环境变量，注意numpy库在/lib64/python2.7/site-packages目录中的

1
export PYTHONPATH=/home/ctan/python2/numpy_installation/lib64/python2.7/site-packages:$PYTHONPATH

ok啦，以上是题外话，安装完numpy之后就可以运行脚本了，将POSCAR和CONTCAR的坐标从分数坐标转换为笛卡尔坐标

POSCAR

1
2
3
4
5
6
7
8
9
10
11
12
13
14
[ctan@baifq-hpc141 primitive-slab-opt-1]$ cat POSCAR_C
CONTCAR\(1\1\1)
1.0
        2.5701000690         0.0000000000         0.0000000000
       -1.2850500345         2.2257719501         0.0000000000
        0.0000000000         0.0000000000        21.2954006195
   Cu
    4
Selective
Cartesian
        +0.0000000000   +0.0000000000   +0.0000000000  F F F
        +1.2850628850   +0.7419165641   +2.0984487770  F F F
        -0.0000128505   +1.4838553860   +4.1968975541  T T T
        +0.0000000000   +0.0000000000   +6.2953463311  T T T

CONTCAR_C

1
2
3
4
5
6
7
8
9
10
11
12
13
14
[ctan@baifq-hpc141 primitive-slab-opt-1]$ cat CONTCAR_C
CONTCAR\(1\1\1)
   1.00000000000000
     2.5701000690000000    0.0000000000000000    0.0000000000000000
    -1.2850500345000000    2.2257719501000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.2954006195000005
  Cu/0e71e558f37
               4
Selective
Cartesian
        +0.0000000000   +0.0000000000   +0.0000000000  F F F
        +1.2850628850   +0.7419165641   +2.0984487770  F F F
        -0.0000084856   +1.4838528659   +4.2004787017  T T T
        -0.0000015980   +0.0000009226   +6.2855788488  T T T

Alright, I’m back. 突然想起昨天球磨的样品今天结束了，连忙赶回实验室取出来装进样品管里等待封管。

By the way，我真是和鱼一样，记忆只有七秒、、、😅。师兄教的球磨罐和球珠的洗涤方法忘记下来了，今天不出意外地给忘了😂。汲取教训，每次出实验室后都更加努力地总结和复盘，仔细回想每一个实验细节

优化前后对比

Alright，扯远了，分析一下顶层和亚顶层优化前后地结构变化吧：

  • POSCAR: (+6.2953463311) - (+4.1968975541) = 2.0984487769999998
  • CONTCAR: (+6.2855788488) - (+4.2004787017) = 2.0851001471000004

前后变化了：(2.0851001471000004-2.0984487769999998)/2.0984487769999998 = -0.006361189296735172

可以看到，表层原子向体相收缩

能量变化

1
2
3
4
5
6
[ctan@baifq-hpc141 primitive-slab-opt-1]$ grep '  without' OUTCAR
  energy  without entropy=      -13.96941183  energy(sigma->0) =      -13.97673168
  energy  without entropy=      -13.96945528  energy(sigma->0) =      -13.97676958
  energy  without entropy=      -13.97006413  energy(sigma->0) =      -13.97726241
  energy  without entropy=      -13.97020643  energy(sigma->0) =      -13.97735554
  energy  without entropy=      -13.97017969  energy(sigma->0) =      -13.97735013

能量从-13.97673168变化为-13.97735013，降低了0.000618449999999271 eV，即弛豫过程是放热，刚切开的表面不稳定，表层原子向体相收缩后，体系能量降低，变得更稳定

小Tips：可以使用python环境便捷地进行基本运算

python环境下进行基本运算

从零开始的科研之旅 > VASP
#化学 #科研
以p(1x1)Cu(111)为例分析表面弛豫的结果
https://hydrogen1222.xyz/2025/04/27/科研/VASP/Starry Sky/以p(1x1)Cu(111)为例分析表面弛豫的结果/
作者
Storm Talia
发布于
2025年4月27日
许可协议