lammps新手水合物分子动力学模拟
我用lammps来模拟PETG的拉伸试验，无论这个程序跑多久，erate值是多少，链始终拉不断。即使在ovito中观察到链已经被拉到足够长，他还是不断。所以我的问题就是怎么才能把它拉断。下面是我的源文件：variable projectindex "PETG" # project name
variable projectindex "petg" # project name
variable sourceindex ../build # data directory
variable paramsindex ../build # parameter directory
variable temperature index 300# system temperature
variable tdampindex 100# temperature damping
variable dielectric index 1# medium dielectric
variable kappaindex 4# electrostatics kappa
variable cutoffindex 9.5# standard cutoff
variable charge_cutoff index 9.5# charge cutoff
variable precision index 0.001# kspace precision
variable lseedindex 723853# langevin seed
variable vseedindex 1486234# velocity init seed
variable tequilindex 50000# equilibration time
variable trunindex 10000000 # run time
variable frestart index 0# 0: equil, 1: restart
variable dtrestart index 100000# delta restart time
variable dtdumpindex 100000# delta dump time
variable dtthermo index 1000# delta thermo time
variable timestep index 0.5# integration time step
variable tfreqindex 10# profile sampling freq
variable nsampleindex 1000# profile conf sampling
variable dtimeequal ${tfreq}*${nsample} # profile dtime
variable restartindex ${params}/${project}.restart
# New Variables for output
variable p1equal "time"
variable p2equal "vol"
variable p3equal "pe"
variable p4equal "temp"
variable p5equal "etotal"
variable p6equal "density"
variable p7equal "mol"
if "${frestart} != 0" then &
"variable dataindex ${restart}" &
else &
"variable dataindex ${project}.data" &
# LAMMPS atomistic input script
echoscreen
unitsreal
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}
fixmom all momentum 100 linear 1 1 1 angular
# Set thermo output
thermo1000
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
computemsd1 all msd/chunk cc1
#fix msd1 all ave/time 1000 1 1000 c_msd1file tmp.out mode vector
#########################################################
# Equilibration #
#########################################################
# STAGE1 NVT from 1000 to 300 K
reset_timestep 0
timestep ${timestep}
#velocityall create ${temperature} 1231
fix1 all nvt temp 1000 300 100 drag 0.1
run5000
unfix 1
# STAGE2 NPT from 300K to 300K
fix 1 all npt temp 300 300 100 iso 0 0 1000 drag 0.1
run5000
unfix1
# Output file
fixdef1 all print 1000 "${p1} ${p2} ${p3} ${p4} ${p5} ${p6} ${p7}" fileequilibration.txt screen no
dump1 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
###############################################
###############################################
variablesimname equal PETG
variable equisteps equal 50000
variable thermosteps equal 100
variable strain_rate equal 1e-5
###############################################
# Deformation
reset_timestep0
thermo ${thermosteps}
timestep 0.5
fix1 all npt temp 300 300 25 y 0.0 0.0 500 z 0.0 0.0 500 nreset 1000
fix2 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"
variablep1 equal "step"
variablep2 equal "dt"
variablep3 equal "pe"
variablep4 equal "ke"
variablep5 equal "etotal"
variablep6 equal "v_strain"
variablep7 equal "-pxx/10000*1.01325"
variablep8 equal "-pyy/10000*1.01325"
variablep9 equal "-pzz/10000*1.01325"
variablep10 equal "lx"
variablep11 equal "ly"
variablep12 equal "lz"
variablep13 equal "press"
variablep14 equal "pxy"
variablep15 equal "pxz"
variablep16 equal "pyz"
variablep17 equal "temp"
variablet2 equal "epair"
variablet3 equal "ebond"
variablet4 equal "eangle"
variablet5 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
run20000
write_restart ${project}.restart2
write_data data.*
print "All Done!"
我问过我的导师，他说可以用一下fix bond/react来试一下，可是我发现那个指令需要拓扑结构（molecule指令）。
（语句表达比较混乱，望见谅） 用的力场就不支持研究化学键的断裂，所以无论你再怎么拉都不会断裂 :handshake
页:
[1]