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的话，模型的体积会大幅度波动，膨胀收缩很剧烈。