球形石墨烯建模教程

背景介绍

石墨烯因为优异的物理和力学性质，是一种很重要的二维材料。其中一个用途就是通过和金属复合来提高基体的力学性能。但是传统的层状复合，存在一些难以避免的问题，比如受载荷的强各向异性，石墨烯本身的团聚等，都会显著降低其实际应用价值。于是，构建三维化的石墨烯网格成为一种潜在的解决手段。包括石墨烯/碳纳米管的结合，蜂窝状石墨烯，网状石墨烯，以及今天我主要关心的球状石墨烯。

球状石墨烯，最简单的例子是足球烯，也就是常说的C60,这一类球形的，完全由近似正六边形和正五边形组成的结构，还有一个名字叫做Goldberg polyhedron (GP)。这类结构有两个重要参数m,n,代表从任意一个正五边形出发，直走m步,左转60°n步，就能刚好到达下一个正五边形，简称G(m,n)。比如下图就是G(1,4).
G(1,4) 图形
这里，构建球形石墨烯，我们就建立这样的初始构型来表征其初始原子分布。

模型建立

一般来说，建立模型我都是用python来做的，尤其是晶体结构比较明确的那种，比如FCC晶体，单层石墨烯这种，不过今天这个由于比较复杂，所以我认为应该是有现成的轮子可以用的。一通谷歌以后，还真找到了几个：

这两个软件都可以用来生成GP,NanoCap自带图形界面，windows/ubuntu都可以直接使用，但是只有一个指定参数，就是C原子数目，它并不是按照数学上的方式来直接建立GP，似乎是加的有优化，好处是生成的球状石墨烯结构上很稳定，坏处就是6000个原子以上我就生成不出来了，也就是规模十分有限。今天重点来讲解Antiprism使用。

Antiprism使用

它的官网显示windows下可以直接使用，但是当你下载.exe文件时候，好家伙，直接没有，无奈我先是安装了他的python包，安装完以后使用过程的语法，好像是有点不对劲，生成的东西不对，也可能是我没有用对吧，最后我是用他的linux版本生成的。真的，这个软件的文档好难理解，甚至我感觉有些像是隐藏参数，so hard.

安装linux版本的Antiprism
1
2
3
sudo add-apt-repository ppa:antiprism/ppa
sudo apt-get update
sudo apt-get install antiprism
安装过程中显示有些包没有更新或者已经没有用不影响最终的使用，我这里是使用WSl来安装的，没有可视化。
生成一个较大的GP(27,25)
1
geodesic -c 27,25 ico | pol_recip | canonical > G27-25.off
其中27，25就是m和n，这几个命令具体意思可以在官网查看，不过真的不是很好理解。如下图这样就是运行成功了。

转换生成的.off文件为MD文件
这里我们只需要他的坐标文件，而且需要进行缩放来使得初始的C-C键长大致在0.142nm，我写了一个简单的python脚本来转换为data文件，这种格式在LAMMPS中很方便使用。

import freud # conda install -c conda-forge freud 很好用的库
import numpy as np
import pandas as pd

def off2data(iname):
    file = open(iname).readlines() # 读取文件
    N = int(file[1].split()[0])
    pos = np.array([i.split() for i in file[2:2+N]],dtype=float) # 获取原始坐标
    Lbox = np.max(pos,axis=0)-np.min(pos, axis=0)
    box = freud.box.Box(*Lbox*1.5)
    nlist = freud.locality.AABBQuery(box, pos).query(pos, dict(num_neighbors=3, exclude_ii=True)).toNeighborList()
    avebond = np.mean(nlist.distances) # 计算平均C-C键长
    factor = 1.42/avebond
    pos *=factor # 缩放原始坐标
    pos -= np.mean(pos, axis=0) # 质心平移至原点
    df = pd.DataFrame(pos, columns=['x', 'y', 'z'])
    df['type'] = '1'
    df['idx'] = np.arange(len(df))+1
    df = df[['idx','type', 'x', 'y', 'z']]
    lo = np.min(pos, axis=0)-20
    hi = np.max(pos, axis=0)+20
    oname = iname.split('.')[0] + '.data'
    with open(oname, 'w') as op: # 写入data文件
        op.write(f'# {oname} written by HerrWu.\n')
        op.write(f'{len(pos)} atoms\n1 atom types\n')
        op.write(f'{lo[0]} {hi[0]} xlo xhi\n')
        op.write(f'{lo[1]} {hi[1]} ylo yhi\n')
        op.write(f'{lo[2]} {hi[2]} zlo zhi\n\n')
        op.write('Atoms # atomic\n\n')

    df.to_csv(oname, sep = ' ',index=False, header=None, mode='a')

if __name__ == '__main__':
    iname = 'G27-25.off'
    off2data(iname)

使用Ovito可视化生成的构型
使用Ovito打开.data文件以后，使用create_bonds，这里对于五边形的部分应该单独做处理的，我就简单直接看了。
全局示意图

初始构型弛豫

使用LAMMPS来进行初始的构型弛豫，in文件如下：

dimension 3
units metal
boundary p p p
atom_style atomic

read_data ./G27-25.data 

neighbor 1.0 bin
neigh_modify every 5 delay 0 check yes

mass 1 12.010700
pair_style airebo/morse 3.0 1 
pair_coeff * * ./CH.airebo-m C

timestep 0.001
shell mkdir minimize
dump mini all custom 100 minimize/mini.*.dump id type x y z vx vy vz
min_style cg
minimize 1.0e-8 1.0e-8 10000 10000
undump mini 
reset_timestep 0 

variable Press equal press
variable Temp equal temp
variable Step equal step
variable PE equal pe

thermo 100
thermo_style custom step temp press pe ke

velocity all create 300 429816734892 dist gaussian mom yes rot yes

fix 1 all nvt temp 300 300 0.5

fix 2 all print 100 "${Step} ${Temp} ${Press} ${PE}" file relax.txt screen no title "Step Temp Press PE"
shell mkdir relax
dump 1 all custom 100 relax/relax.*.dump id type x y z vx vy vz

run 20000