博文

LAMMPS中使用fix deform命令拉伸单层石墨烯的算例子+附in文件

已有 21242 次阅读 2018-7-23 20:00 |系统分类:科研笔记

fix defom +fix npt 来实现单轴拉伸单层石墨烯 
前几天把很久以前写的单轴拉伸的in文件拿出来整理了一下，发现in文件的可读性并不强
所以就重新编写了一下，计算的目的是单轴拉伸石墨烯片，in文件如下:

# 1------general---------------
units metal      #单位为lammps 中的metel 类型
dimension3          #模拟的维度 三维
boundaryp p p      #周期边界条件 
atom_styleatomic   #原子类型自动
atom_modify map array sort 100 2.0
neighbor 3.0 bin
neigh_modify every 1 delay 0 check yes
#2 -----global variable-------定义全局变量，方便统一修改；
variable temperature equal 0.1
variable tstep equal 0.001
variable pressure equal 0
variable thermalstep equal 100
variable dumpstep equal 1000
variable relaxtime equal 100000
variable totaltime equal 500000
variable deformrate equal 10e-4
#3 -------model------------
#box tilt large
read_data      graphene.data
#4------potential----------
pair_style     airebo 3.0 1 1
pair_coeff     * * CH.airebo C  
# 5------minimize -----------
min_style      sd
minimize       1.0e-4 1.0e-6 100 1000
#6----define computes------
variable nums equal count(all) #计算原子数目
variable toval equal ly*lz*3.35 #计算系统总体积，单层石墨烯厚度为0.335nm
variable vol equal ${toval} #体系的初始体积
variable vatom equal v_vol/v_nums #单个原子体积
compute 1 all stress/atom NULL #计算体系中单原子应力
compute xalls all reduce sum c_1[1]
variable xstress equal (c_xalls)/(v_toval*10000)
compute yalls all reduce sum c_1[2]
variable ystress equal (c_yalls)/(v_toval*10000)
compute zalls all reduce sum c_1[3]
variable zstress equal (c_zalls)/(v_toval*10000)
#7---------relaxation-------------------------
timestep ${tstep}
velocity all create ${temperature} 231743 units box dist gaussian
dump   1 all atom 100 all-cy.lammpstrj
fix    1 all npt temp 0.1 0.1 0.1 aniso 0 0 10 drag 1  #三个方向独立控压
compute  2 all pe/atom     
compute  pe all reduce sum c_2
variable PE equal "c_pe"
fix    pe all print 100 "${PE} " file PE.txt screen no
thermo  ${thermalstep}
#thermo_style  custom step pe ke etotal lx ly lz
run   ${relaxtime}
#8 -------tension-----------------------------
reset_timestep 0
unfix  1
undump  1
unfix  pe
fix avestress all ave/atom 1 ${dumpstep} ${dumpstep} c_1[1] c_1[2] c_1[3] c_1[4] c_1[5] c_1[6]

#-----calculation of strain and Cumulative compression--------
variable                tmp equal "ly"
variable                L0 equal ${tmp}
variable                strain equal "(ly - v_L0)/v_L0"
variable                Cumulativels equal "(ly - v_L0)/10"  
fix Step all print 100 "${strain} ${Cumulativels} ${ystress}" file grapoten.txt screen no

dump 2 all custom ${dumpstep} tension.xyz type x y z f_avestress[1] f_avestress[2] 
f_avestress[3] f_avestress[4] f_avestress[5] f_avestress[6]

fix  1 all npt temp ${temperature} ${temperature} 0.1 x ${pressure} ${pressure} 10 
z ${pressure} ${pressure} 10

thermo_style custom step pe ke etotal lx ly lz
fix   2 all deform 1 y erate ${deformrate} remap x units box 
fix   3 all ave/time 2 500 1000 v_xstress v_ystress v_zstress file pressure.out
run   ${totaltime}

这次计算出现的主要的问题：在弛豫时候使用fix npt命令对烯片的x、y、z三个方向进行单独的控压

在使用fix npt控压时我们需要指定pdamp参数，这个参数是一个时间单位，pdamp数值不能是任意的太小会导致体系压强和体积的大幅度波动；如果pdamp过大，那么需要经过较长的一段时间压力才能达到平衡。手册中说明对于大多数的模型，pdamp的数值可以取1000倍的时间步长（timestep),本例子中如果按照手册提示应该选择0.001*1000=1，但是经过测试发现如果取1的话，模型的体积会大幅度波动，膨胀收缩很剧烈。