构建用于分子动力学（MD）或者密度泛函（DFT）模拟的氧化石墨烯（GO）模型是一个常见但需要细致处理的任务。通常可以借助Materials Studio的可视化界面帮助建模。但是如果构建的尺寸比较大的话，MS手动建模也显得有点低效。我们在这里编写了一个python小程序给大家参考如何构建大尺寸氧化石墨烯片模型。
碳基底上共价键含氧官能团，主要包括：
- 环氧基：位于碳环上方或下方。
- 羟基：与碳原子直接键合。
- 羧基：位于石墨烯片的边缘。本次我们介绍周期性结构的氧化石墨烯片，只考虑羟基的情况。
Python脚本的基本逻辑是：确定石墨烯片的大小尺寸，构建初始石墨烯结构。然后随机选择一定数量的碳原子。对于羟基：选择一个碳原子，然后连接一个OH官能团。 运行完脚本之后，得到结构文件OUT.vasp。此格式文件可以用VESTA可视化软件打开(www.jp-minerals.org/vesta/en/ )，并可另存为其它格式文件（如cif等）。
石墨烯的晶格常数为2.46埃，原胞为六角格子，但是在构建石墨烯片的时候，可以选择更大的矩形超胞。
python PyGO.py 7 12 25 745236
# 7 和 12 分别为x和y方向的基矢倍数
# 25 为需要添加的OH官能团的数量
# 745236 为随机数种子
运行完之后，用vesta打开OUT.vasp文件，将会看到构建好的氧化石墨烯片，如下图。
# By 南方小土豆 https://juejin.cn/post/7542769933726433295
import numpy as np
import os, sys
inputs = sys.argv
bb = 2.46
aa = np.sqrt(3.0)*bb
h = 10
def get_dist(AA,BB,la,lb):
ss = 10
for pp in BB:
for ii in range(-1,2):
lx = ii*la
for jj in range(-1,2):
ly = jj*lb
tt = (AA[0]-pp[0]+lx)**2 + (AA[1]-pp[1]+ly)**2 + (AA[2]-pp[2])**2
if np.sqrt(tt) < ss:
ss = np.sqrt(tt)
return ss
def get_dd(AA,pp):
ss = (AA[0]-pp[0])**2 + (AA[1]-pp[1])**2 + (AA[2]-pp[2])**2
return np.sqrt(ss)
def get_two_posc(AA,BB,la,lb):
nd = []
ss = 1.6
ids = 0
xx = []
yy = []
for pp in BB:
ids+=1
for ii in range(-1,2):
lx = ii*la
for jj in range(-1,2):
ly = jj*lb
tt = (AA[0]-pp[0]+lx)**2 + (AA[1]-pp[1]+ly)**2 + (AA[2]-pp[2])**2
if np.sqrt(tt) < ss and np.sqrt(tt)>1.0:
nd.append(ids-1)
xx.append((AA[0]+pp[0]-lx)/2.0)
yy.append((AA[1]+pp[1]-ly)/2.0)
return nd[0],xx[0],yy[0]
def get_gra(nx,ny,noh,seeds):
np.random.seed(seeds)
p1 = [0,0,0.5*h]
p2 = [1./3.*aa,0,0.5*h]
p3 = [0.5*aa,0.5*bb,0.5*h]
p4 = [5.0/6*aa,0.5*bb,0.5*h]
pos0 = np.array([p1,p2,p3,p4])
pos_c =[]
for ii in range(nx):
la = ii*aa
for jj in range(ny):
lb = jj*bb
for pp in pos0:
pos_c.append([pp[0]+la, pp[1]+lb, pp[2]])
c_label = []
######## OH
pos_oh_o=[]
pos_oh_h=[]
ncc=0
ss = 0
for ii in range(20*noh):
idc = int(len(pos_c)*np.random.rand())
if idc not in c_label and ncc<noh:
xr = np.random.rand()
zz = 0.5*h+1.5
if xr<0.5:
zz = 0.5*h-1.5
if len(pos_oh_o)>0:
ss = get_dist([pos_c[idc][0], pos_c[idc][1],zz ], pos_oh_o, nx*aa, ny*bb )
if (len(pos_oh_o)>0 and ss>2.5) or len(pos_oh_o)==0:
pos_oh_o.append([pos_c[idc][0], pos_c[idc][1],zz ])
pos_oh_h.append([pos_c[idc][0], pos_c[idc][1],zz+np.sign(xr-0.5)])
c_label.append(idc)
ncc+=1
ele = ['C','O','H']
num = [len(pos_c),len(pos_oh_o), len(pos_oh_h)]
f = open('OUT.vasp','w')
print('#POSCAR graphene oxide', file=f)
print('1.0 ', file=f)
print( '%12.6f %12.6f %12.6f' %(nx*aa, 0, 0), file=f )
print( '%12.6f %12.6f %12.6f' %(0, ny*bb, 0), file=f )
print( '%12.6f %12.6f %12.6f' %(0, 0, h), file=f )
for ii in range( len(ele) ):
key = ele[ii]
if ii==len(ele)-1:
print(key, file=f)
else:
print(key, end=' ', file=f)
for ii in range( len(num) ):
key = num[ii]
if ii==len(num)-1:
print(key, file=f)
else:
print(key, end=' ', file=f)
print('Cartesian ', file=f)
for pp in pos_c:
print( '%12.6f %12.6f %12.6f' %(pp[0], pp[1], pp[2]), file=f)
for pp in pos_oh_o:
print( '%12.6f %12.6f %12.6f' %(pp[0], pp[1], pp[2]), file=f)
for pp in pos_oh_h:
print( '%12.6f %12.6f %12.6f' %(pp[0], pp[1], pp[2]), file=f)
f.close()
return
nx=int(inputs[-4])
ny=int(inputs[-3])
noh=int(inputs[-2])
seeds=int(inputs[-1])
get_gra(nx,ny,noh,seeds)
print('Lx= %6.2f, Ly= %6.2f' %(nx*aa,ny*bb) )
