24小时热门版块排行榜

返回列表

dongzwang

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 3.5
帖子: 7
在线: 1.5小时
虫号: 7773077
注册: 2018-01-15

[交流] 做高分子模拟拉伸试验，链始终不断已有3人参与

我用lammps来模拟PETG的拉伸试验，无论这个程序跑多久，erate值是多少，链始终拉不断。即使在ovito中观察到链已经被拉到足够长，他还是不断。所以我的问题就是怎么才能把它拉断。下面是我的源文件：
variable project  index "PETG" # project name
variable project  index "petg" # project name
variable source  index ../build # data directory
variable params  index ../build # parameter directory
variable temperature index 300  # system temperature
variable tdamp  index 100  # temperature damping
variable dielectric index 1  # medium dielectric
variable kappa  index 4  # electrostatics kappa
variable cutoff  index 9.5  # standard cutoff
variable charge_cutoff index 9.5  # charge cutoff
variable precision index 0.001  # kspace precision
variable lseed  index 723853  # langevin seed
variable vseed  index 1486234  # velocity init seed
variable tequil  index 50000  # equilibration time
variable trun  index 10000000 # run time
variable frestart index 0  # 0: equil, 1: restart
variable dtrestart index 100000  # delta restart time
variable dtdump  index 100000  # delta dump time
variable dtthermo index 1000  # delta thermo time
variable timestep index 0.5  # integration time step
variable tfreq  index 10  # profile sampling freq
variable nsample  index 1000  # profile conf sampling
variable dtime  equal ${tfreq}*${nsample} # profile dtime
variable restart  index ${params}/${project}.restart
# New Variables for output
variable p1  equal "time"
variable p2  equal "vol"
variable p3  equal "pe"
variable p4  equal "temp"
variable p5  equal "etotal"
variable p6  equal "density"
variable p7  equal "mol"
if "${frestart} != 0" then &
"variable data  index ${restart}" &
else &
"variable data  index ${project}.data" &

# LAMMPS atomistic input script

echo  screen
units  real
atom_style full

# Interaction potential definition

pair_style lj/class2/coul/long ${cutoff} ${charge_cutoff}
bond_style harmonic
special_bonds lj/coul 0 0 1
if "${frestart} != 0" then "read_restart ${data}" else "read_data ${data}"
include  ${project}.params

# Integration conditions (check)

timestep ${timestep}
kspace_style pppm/cg ${precision}
dielectric ${dielectric}
fix  mom all momentum 100 linear 1 1 1 angular

# Set thermo output
thermo  1000
thermo_style custom step pe etotal temp time vol density

# Calculate the mean squared displacement for Multiple chunks of atoms
compute    cc1 all chunk/atom molecule
compute  msd1 all msd/chunk cc1
#fix msd1 all ave/time 1000 1 1000 c_msd1

file tmp.out mode vector

#########################################################
# Equilibration                #
#########################################################
# STAGE1 NVT from 1000 to 300 K
reset_timestep 0
timestep ${timestep}
#velocity  all create ${temperature} 1231
fix  1 all nvt temp 1000 300 100 drag 0.1
run  5000
unfix 1
# STAGE2 NPT from 300K to 300K
fix       1 all npt temp 300 300 100 iso 0 0 1000 drag 0.1
run  5000
unfix  1

# Output file
fix  def1 all print 1000 "${p1} ${p2} ${p3} ${p4} ${p5} ${p6} ${p7}" file  equilibration.txt screen no
dump  1 all cfg ${dtdump} equi.dump.* mass type xs ys zs id mol q x y z mol
#dump_modify    1 element o_2 ho2 oh c hc c1 c_1 o_1 cp c_0 ho

###############################################
###############################################

variable  simname equal PETG
variable equisteps equal 50000

variable thermosteps equal 100
variable strain_rate equal 1e-5

###############################################
# Deformation

reset_timestep  0
thermo       ${thermosteps}
timestep 0.5
fix  1 all npt temp 300 300 25 y 0.0 0.0 500 z 0.0 0.0 500 nreset 1000
fix  2 all deform 1 x erate ${strain_rate} units box remap x

variable tmp equal "lx"
variable L0 equal ${tmp}
variable strain equal "(lx - v_L0)/v_L0"
variable  p1 equal "step"
variable  p2 equal "dt"
variable  p3 equal "pe"
variable  p4 equal "ke"
variable  p5 equal "etotal"
variable  p6 equal "v_strain"
variable  p7 equal "-pxx/10000*1.01325"
variable  p8 equal "-pyy/10000*1.01325"
variable  p9 equal "-pzz/10000*1.01325"
variable  p10 equal "lx"
variable  p11 equal "ly"
variable  p12 equal "lz"
variable  p13 equal "press"
variable  p14 equal "pxy"
variable  p15 equal "pxz"
variable  p16 equal "pyz"
variable  p17 equal "temp"
variable  t2 equal "epair"
variable  t3 equal "ebond"
variable  t4 equal "eangle"
variable  t5 equal "edihed"

fix def1 all print 100 "${p1} ${p6} ${p7} ${p8} ${p9} ${p10} ${p11} ${p12} ${p13} ${p14} ${p15} ${p16} ${p17}" file PETG_npt.def1.txt screen no title "step v_strain pxx pyy pzz lx ly lz press pxy pxz pyz temp"
fix def2 all print 100 "${p1} ${t2} ${t3} ${t4} ${t5}" file PETG_npt.def2.txt screen no
fix myprint all print 100 "${p1} ${p2} ${p3} ${p4} ${p5}" file PETG_npt.engergy.txt screen no title "step dt pe ke etotal"
dump 2 all custom 10000 dump_npt.mole.* id type mol fx fy fy mass x y z

thermo_style custom step dt pe ke etotal pxx pyy pzz lx ly lz epair ebond eangle edihed

run  20000
write_restart ${project}.restart2
write_data data.*

print "All Done!"
我问过我的导师，他说可以用一下fix bond/react来试一下，可是我发现那个指令需要拓扑结构（molecule指令）。
（语句表达比较混乱，望见谅）

回复此楼