Hello, LAMMPS community!
I have a reason to think that LAMMPS fix bond/break either does not break bonds correctly or does not report the breakage correctly.
I perform 1D stretching of a simulation box filled with a Kremer-Grest polymer network with breakable bonds. The bonds are of two types : 1 and 2 - mimicing polymer backbone bonds and sulfur cross-link bonds.
LAMMPS command in the terminal is [ do not ask, please, why I do not use the LAMMPS job submission script with srun routine and #SBATCH directives ] :
mpirun -np 112 <LAMMPS executable> -in kg_1D_str_X-lnk_BB_brk.lam -var input <LAMMPS input file> -var base <base name> -var tRlx 500 -log <log file>
LAMMPS configuration script kg_1D_str_X-lnk_BB_brk.lam looks like as follows :
# Uniaxial stretching a KG model network from ElF 1 to 20 along X axis
# X-links & backbone bonds are breakable with the distance criteria
# #################
# ### VARIABLES ###
# #################
variable mass equal 1.0 # m_b -- mass of KG bead
variable tDamp equal 2.0 # tau -- relaxation time of temperature for Langevin dynamics
variable T equal 1.0 # final temperature for Langevin dynamics
variable tStep equal 0.01 # time step
variable finElF equal 20 # target elongation factor (ElF)
variable tElF equal ${tRlx} # time per strain increment
variable strRate equal 1/${tElF} # strain rate (engineering)
variable nRun equal &
(${finElF}-1)*${tElF}/${tStep} # steps to make ElF equal to finElF
# print deformation mode info
print " Final elongation factor = ${finElF}"
print " Time per strain increment = ${tElF} tau"
print " Strain rate = ${strRate} tau^-1 (engineering)"
print " Number of time steps = ${nRun}"
# ######################
# ### INITIALISATION ###
# ######################
units lj # LJ units
dimension 3 # 3D problem dimension
boundary p p p # periodic BCs
atom_style molecular # molecular atom style
# ######################################
# ### NEIGHBOUR LIST & COMMUNICATION ###
# ######################################
neighbor 2.0 bin # neighbour skin thickness & type
neigh_modify every 1 delay 0 check yes # neighbour list modification settings
comm_modify mode single cutoff 4.0 # inter-process commmunication
# ##################
# ### READ INPUT ###
# ##################
read_data ${input}
mass * ${mass} # KG bead mass
# ####################
# ### INTERACTIONS ###
# ####################
variable cutoff equal 1.12246204830937
pair_style lj/cut ${cutoff}
pair_coeff * * 1.0 1.0 ${cutoff}
pair_modify shift yes
bond_style fene
bond_coeff * 30.0 1.5 0.0 1.0
special_bonds lj 1 1 1
angle_style cosine
angle_coeff * 0.2060
# ################
# ### DYNAMICS ###
# ################
fix dynamics all langevin ${T} ${T} ${tDamp} 432985345 gjf vfull tally yes
fix integrator all nve
# ##############################
# ### STRETCH BOX UNIAXIALLY ###
# ### PRESERVING VOLUME ###
# ##############################
fix 1DStretching all deform 10 &
x erate ${strRate} &
y volume &
z volume &
units box &
remap x
# rebalance grid of nodes when box is strained
fix rebal all balance 10000 0.95 shift xyz 10 1.05
# ######################################
# ### BREAK BACKBONE BONDS & X-LINKS ###
# ######################################
fix brkBB all bond/break 10 1 1.3
fix brkX all bond/break 10 1 1.275
# #############
# ### FIXES ###
# #############
variable lambda equal lx/vol^(1.0/3.0) # x-ElF
# number of broken bonds
fix brkAvg all ave/time &
1000 1 1000 v_lambda f_brkBB[1] f_brkX[1] f_brkBB[2] f_brkX[2] &
file brkNum.t
# ######################
# ### TIMESTEP & RUN ###
# ######################
timestep ${tStep}
run ${nRun} upto start 0 stop ${nRun} every 50000 NULL
As you can see, there is an accidentally made typo, namely :
fix brkBB all bond/break 10 1 1.3
fix brkX all bond/break 10 1 1.275
i.e. breakage is defined only for the bonds of type 1, not for the type 2. My bad, I will fix it for a new MD run. Nonetheless, I have noticed that brkNum.t outpuf file has the following interesting data :
# Time-averaged data for fix brkAvg
# TimeStep v_lambda f_brkBB[1] f_brkX[1] f_brkBB[2] f_brkX[2]
0 1 0 0 0 0
1000 1.02 0 0 0 1
2000 1.04 0 0 0 1
3000 1.06 0 0 0 1
...
243000 5.86 0 1 660 7294
244000 5.88 0 0 670 7400
245000 5.9 0 0 679 7497
...
i.e. the type 2 bonds are also broken even though for them breakage is not specified. Thus, there are two questions I would like to address :
-
how come LAMMPS
fix bond/breakactually breaks type 2 bonds? -
how come the cumulative number of broken bonds got increased, i.e.
1000 1.02 0 0 0 1, whereas there was not a single bond broken on any breakage time step?
P.S. Does any body know how to output info about broken bonds, e.g. their coordinates or the coordinates of beads that a broken bond connect? Thank you!
I’m sorry for making such a silly mistake. Thanks for questioning, this helped to find the truth as always.