AIREBO势函数截断半径测试

背景介绍

AIREBO势函数是常用的模拟石墨烯，碳纳米管等结构的势函数，其近程相互作用默认在0.17nm-0.2nm之间会有一个非物理的应变硬化，已经有一些研究手动修改截断半径到0.2nm，比如这篇.这里我们简单记录一下测试不同的截断半径对石墨烯单轴拉伸曲线的影响。

建立单层石墨烯

这里为了计算快速就建立一个5nmx5nm的小石墨烯，bulid_structure函数见文章.

def write_data(pos, La, Lb, Lc, filename='grap.data'):
    print(f'Writing {filename}...')
    xmin, ymin, zmin = np.min(pos, axis=0)
    
    with open(filename, 'w') as op:
        # First line is a comment line
        op.write('Graphene structure, written by HerrWu\n\n')
        op.write('{} atoms\n'.format(len(pos)))
        op.write('1 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-Lc, zmin+Lc))

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

        # Write each Position
        aidx = 1
        for i in pos:
            op.write('{} 1 {} {} {}\n'.format(aidx, *i))
            aidx += 1
            
x,y,z = 12, 21, 1
Lattice_Constant = 1.42
Structure_Type = 'GRA'
pos = bulid_structure(Lattice_Constant, Structure_Type, x, y, z)
La, Lb, Lc = x*Lattice_Constant*3, y*Lattice_Constant*3**0.5, 10
write_data(pos, La, Lb, Lc, filename='grap.data')

石墨烯拉伸脚本

这个脚本是直接从网上抄的。

#uniaxial tensile test of graphene
##---------------INITIALIZATION-------------------------------
units          metal
dimension 	    3 
boundary       p p f
atom_style 	atomic
newton 		on
##---------------ATOM DEFINITION------------------------------
read_data 	grap.data
##---------------FORCE FIELDS---------------------------------
pair_style 	airebo 3.0
pair_coeff     * * CH.airebo C
##---------------SETTINGS-------------------------------------
timestep 	0.0005
variable   ts equal 0.0005
##---------------COMPUTES-------------------------------------
compute 	1 all stress/atom NULL
compute    2 all reduce sum c_1[1] c_1[2]

variable   Lx equal lx
variable   Ly equal ly
variable   Lz equal lz
variable   Vol equal vol
variable   thickn equal 3.4
fix		1 all npt temp 300 300 0.05 x 0 0 0.5 y 0 0 0.5
thermo 	2000
##---------------RELAXATION--------------------------------------
run            50000
##---------------DEFORMATION-------------------------------------
unfix              1
reset_timestep     0
fix		1 all npt temp 300 300 0.05 x 0 0 0.5
fix        2 all ave/time 1 100 100 c_2[1] c_2[2]
fix        3 all ave/time 1 100 100 v_Lx v_Ly v_Lz v_Vol
variable   srate equal 1.0e9
variable   srate1 equal "v_srate / 1.0e12"
fix		4 all deform 1 y erate ${srate1} units box remap x
run            100
##---------------THERMO-OUTPUTS-----------------------------------
variable   CorVol equal f_3[4]*v_thickn/(f_3[3])
variable   ConvoFac equal 1/1.0e4
variable   sigmaxx equal f_2[1]*v_ConvoFac/v_CorVol
variable   sigmayy equal f_2[2]*v_ConvoFac/v_CorVol
variable   StrainPerTs equal v_srate1*v_ts
variable   strain equal v_StrainPerTs*step
thermo 	100
thermo_style custom step temp v_strain v_sigmaxx v_sigmayy pe ke lx ly vol 
##---------------DEFORMATION-------------------------------------
dump           1 all atom 5000 tensile_test.lammpstrj
run            500000