石墨烯网络与铝基体复合建模教程

背景介绍

上一篇文章中我们建立了一个球状的石墨烯，这里我们对其进行一定的空间排布，然后与金属铝混合在一起，为之后计算其力学性质做准备。

BCC空间排布

BCC是一种常见的金属晶体结构，比如Fe，当然，这里的排布方式并不是固定了，也可以采用FCC或者简单的堆垛，采用BCC我目前也没有必须的科学依据。我写的python代码如下。

import pandas as pd
import numpy as np

def getBCCgra(iname):
    data = pd.read_csv('G27-25.data', skiprows=9, 
                header=None,index_col=False,sep=' ',
                names=['id', 'type', 'x', 'y', 'z'])

    pos = data[['x', 'y', 'z']].values
    pos -= np.mean(pos, axis=0)
    R = np.mean(np.max(pos, axis=0) - np.min(pos, axis=0))/2

    atype = 'BCC'
    Lattice_Constant = (2*R+1.42)*2/np.sqrt(3)
    basis = np.array([[1.0, 0.0, 0.0],
                    [0.0, 1.0, 0.0],
                    [0.0, 0.0, 1.0]]) * Lattice_Constant
    base_atoms = np.array([[0.0, 0.0, 0.0],
                           [0.5, 0.5, 0.5]]) * Lattice_Constant

    x, y, z = 4, 2, 2
    positions = []
    for i in range(x):
        for j in range(y):
            for k in range(z):
                base_position = np.array([i, j, k], dtype=float)
                cart_position = np.dot(basis, base_position)
                for atom_index in range(len(base_atoms)):
                    npos = pos.copy()
                    npos += cart_position + base_atoms[atom_index]
                    for m in npos:
                        positions.append(m)
    positions = np.array(positions)
    xhi,yhi,zhi = Lattice_Constant*(np.array([x,y,z])-1)
    positions = positions[(positions[:,0]<xhi) 
                          & (positions[:,1]<yhi) 
                          & (positions[:,2]<zhi)]
    positions = positions[(positions[:,0]>0) 
                          & (positions[:,1]>0) 
                          & (positions[:,2]>0)]

    df = pd.DataFrame(positions, columns=['x', 'y', 'z'])
    df['type'] = '1'
    df['idx'] = np.arange(len(df))+1
    df = df[['idx','type', 'x', 'y', 'z']]
    oname = iname.split('.')[0] + f'-{x-1}-{y-1}-{z-1}.data'
    with open(oname, 'w') as op: # 写入data文件
        op.write(f'# {oname} written by HerrWu.\n')
        op.write(f'{len(df)} atoms\n1 atom types\n')
        op.write(f'0 {xhi} xlo xhi\n')
        op.write(f'0 {yhi} ylo yhi\n')
        op.write(f'0 {zhi} 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.data'
    getBCCgra(iname)

运行结果如下：

BCC石墨烯弛豫

LAMMPS弛豫脚本如下：

dimension 3
units metal
boundary p p p
atom_style atomic

read_data ./G27-25-3-1-1.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 npt temp 300 300 0.1 x 0 0 1 y 0 0 1 z 0 0 1

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 50000

弛豫后构型球状石墨烯逐渐转变为二十面体，面与面的接触面积增大，整体体积收缩，结果如下：
弛豫后BCC石墨烯网络示意图

石墨烯网络与铝基体复合

首先在盒子内添加FCC铝原子，晶格常数这里取4.05，脚本如下：

import numpy as np
import pandas as pd

def bulid_structure(Lattice_Constant, Structure_Type, x, y, z):
    """
    该函数用于创建常见金属晶体结构！
    :param Lattice_Constant: 晶格常数
    :param Structure_Type: 金属晶体结构：FCC, BCC, HCP
    :param x: x方向拓展个数
    :param y: y方向拓展个数
    :param z: z方向拓展个数
    :return: 原子坐标 N-3 array
    """
    if Structure_Type == 'FCC':
        # Cubic FCC basis
        basis = np.array([[1.0, 0.0, 0.0],
                        [0.0, 1.0, 0.0],
                        [0.0, 0.0, 1.0]]) * Lattice_Constant

        base_atoms = np.array([[0.0, 0.0, 0.0],
                            [0.5, 0.5, 0.0],
                            [0.5, 0.0, 0.5],
                            [0.0, 0.5, 0.5]]) * Lattice_Constant
    elif Structure_Type == 'BCC':
        # Cubic BCC basis
        basis = np.array([[1.0, 0.0, 0.0],
                        [0.0, 1.0, 0.0],
                        [0.0, 0.0, 1.0]]) * Lattice_Constant

        base_atoms = np.array([[0.0, 0.0, 0.0],
                            [0.5, 0.5, 0.5]]) * Lattice_Constant
    elif Structure_Type == 'HCP':
        # Cubic HCP basis
        basis = np.array([[1.0, 0.0, 0.0],
                        [0.0, np.sqrt(3), 0.0],
                        [0.0, 0.0, np.sqrt(8 / 3)]]) * Lattice_Constant

        base_atoms = np.array([[0.0, 0.0, 0.0],
                            [0.5, 0.5 * np.sqrt(3), 0.0],
                            [0.5, np.sqrt(3) * 5 / 6, 0.5 * np.sqrt(8 / 3)],
                            [0.0, 1 / 3 * np.sqrt(3), 0.5 * np.sqrt(8 / 3)]]) * Lattice_Constant
    # Generate atom positions
    positions = []
    for i in range(x):
        for j in range(y):
            for k in range(z):
                base_position = np.array([i, j, k])
                cart_position = np.inner(basis.T, base_position)
                for atom in base_atoms:
                    positions.append(cart_position + atom)
    return np.array(positions)

def build_Al(gra, LAl=4.05):
    xmin, ymin, zmin = np.min(gra, axis=0)
    xmax, ymax, zmax = np.max(gra, axis=0)
    a, b ,c = np.round((np.max(gra, axis = 0) - np.min(gra, axis=0))/LAl).astype(int)
    print(f'Building Al...')
    Al = np.array(bulid_structure(LAl, 'FCC', a, b, c))
    Al = Al - np.min(Al, axis=0) + np.min(gra,axis=0)
    La, Lb, Lc = a*LAl, b*LAl, c*LAl
    atype = np.r_[np.ones(len(Al)), np.ones(len(gra))*2].astype(int)
    return gra, Al, atype, La, Lb, Lc

def write_data(points_mask, atype_mask, La, Lb, Lc, filename='gra-Al.data'):
    print(f'Writing {filename}...')
    xmin, ymin, zmin = np.min(points_mask, axis=0)
    
    with open(filename, 'w') as op:
        # First line is a comment line
        op.write('Gra-Al structure, written by HerrWu\n\n')
        op.write('{} atoms\n'.format(len(points_mask)))
        op.write('2 atom types\n')

        # Specify box Dimensions
        op.write('{} {} xlo xhi\n'.format(xmin, xmin+La))
        op.write('{} {} ylo yhi\n'.format(ymin, ymin+Lb))
        op.write('{} {} zlo zhi\n\n'.format(zmin, zmin+Lc))

        # Atoms section
        op.write('Atoms\n\n')

        # Write each Position
        aidx = 1
        for itype, pos in zip(atype_mask, points_mask):
            op.write('{} {} {} {} {}\n'.format(aidx, itype, *pos))
            aidx += 1

data = pd.read_csv('relax.50000.dump', skiprows=9, 
            header=None,index_col=False,sep=' ',
            names=['id', 'type', 'x', 'y', 'z', 'vx', 'vy', 'vz'])
gra = data[['x', 'y', 'z']].values
gra, Al, atype, La, Lb, Lc = build_Al(gra, LAl=4.05)
points = np.r_[Al, gra]
write_data(points, atype, La, Lb, Lc)

但是此时的铝和石墨烯之间原子存在重叠的现象，需要将距离过近的原子删掉，使用LAMMPS脚本如下：

dimension 3
units metal
boundary p p p
atom_style atomic

read_data gra-al.data 

mass 1 26.981538
mass 2 12.010700

group al type 1
group gra type 2

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes

pair_style hybrid airebo/morse 3.0 1 eam/alloy lj/cut 7.83125 #（2.5*sigma）
pair_coeff * * airebo/morse ../CH.airebo-m NULL C
pair_coeff * * eam/alloy ./Al_DFT.eam.alloy Al NULL
pair_coeff 1 2 lj/cut 0.003457 3.1325

timestep 0.001

delete_atoms overlap 3.1325 al gra # 这句话很重要，再次体现LAMMPS 代码的高效性能 

write_data gra-al-delete.data

删除后模型如下：

最后进行弛豫，LAMMPS脚本如下：

dimension 3
units metal
boundary p p p
atom_style atomic

read_data gra-al-delete.data 

mass 1 26.981538
mass 2 12.010700

group al type 1
group gra type 2

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes

pair_style hybrid airebo/morse 3.0 1 eam/alloy lj/cut 7.83125 #（2.5*sigma）
pair_coeff * * airebo/morse ../CH.airebo-m NULL C
pair_coeff * * eam/alloy ./Al_DFT.eam.alloy Al NULL
pair_coeff 1 2 lj/cut 0.003457 3.1325

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 npt temp 300 300 0.1 x 0 0 1 y 0 0 1 z 0 0 1

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 50000