Hello all,
I am using the 19 Nov 2024 - Development build of LAMMPS and am working on a project attempting to form crosslinking networks of bottlebrush polymers that bond at the endpoints of the side chains. Because my system is relatively small, I want to exclude intramolecular bonds, but when I tried using the molecule keyword in the fix bond/create command (molecule inter), an illegal fix bond/create command error was thrown. Without using this keyword, the simulation runs and produces the default of both inter and intra bonds, but even when using molecule off, the same illegal fix error is thrown. My script is as follows:
--------------- Initialization & System Setup ---------------
units lj
dimension 3
boundary p p p
atom_style full
# Neighbor settings
neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes
comm_modify cutoff 16.0
# ---------------------- Variables (static params) ----------------------
variable ng index 2
variable oh_offset index 0.5
variable backbone_coeff index 0.7
variable sss_coeff index 0.58
variable eps_lo index 0.273 # keep as index so you can override at launch
variable T_hi index 1.0
variable seed index 12345
variable run index 1
# ---------------- Define Atoms, Bonds, Angles (Polymer Topology) ----------------
pair_style lj/cut/coul/long 10.0
bond_style harmonic
special_bonds lj/coul 0.0 0.0 1.0
angle_style harmonic
dihedral_style harmonic
read_data di_retyped_ld_r1_2_80_0.5.data extra/atom/types 1 extra/bond/types 1 extra/bond/per/atom 1 extra/special/per/atom 50
mass * 1.0
group backbone type 1
kspace_style pppm 1e-4
dielectric ${eps_lo}
pair_modify shift yes
pair_coeff * * 1.0 1.0 1.12246
bond_coeff 4 100.0 0.965
# ---------------- Output Settings (Diagnostics and Trajectory) ----------------
compute myRg all gyration
thermo 1000
thermo_style custom step temp pe etotal ecoul c_myRg spcpu
thermo_modify flush yes
dump traj all custom 100000 xlink${run}_ng${ng}_oh${oh_offset}_bb${backbone_coeff}_sss${sss_coeff}_di${eps_lo}_t${T_hi}.lammpstrj id mol type x y z
dump bbtraj backbone custom 100000 xlink${run}_ng${ng}_oh${oh_offset}_bb${backbone_coeff}_sss${sss_coeff}_di${eps_lo}_t${T_hi}.bbtraj id x y z
dump_modify traj sort id
dump_modify bbtraj sort id
log xlink${run}_ng${ng}_oh${oh_offset}_bb${backbone_coeff}_sss${sss_coeff}_di${eps_lo}_t${T_hi}.log
# ---------------- Thermal Equilibration and Pre-anneal Stages ----------------
minimize 1e-4 1e-6 1000 10000
velocity all create {T_hi} 12345 mom yes rot yes dist gaussian
timestep 0.0005
fix 1 all nve/limit 0.05
fix 2 all langevin {T_hi} ${T_hi} 10.0 12345
restart 1000000 restart_xlink_r${run}.*
run 10000
unfix 1
timestep 0.001
fix 1 all nve
fix xlink all bond/create 100 6 6 1.00 4 iparam 1 7 jparam 1 7 prob 0.1 12345 molecule off
run 15000000
write_data xlinked_r1_{ng}\_80\_{oh_offset}.data
and my data file has 6 atom types and 3 bond types, with my script intending to create bonds only between the type 6 atoms, change them to type 7 atoms that don’t bond with other atom types, and classify the new bonds that form as type 4 bonds (for ease of analysis). Can anyone tell at a glance if there is something wrong with my implementation of the fix bond/create command, or whether the molecule keyword is only available to later builds of LAMMPS?
Thanks,
Yuvan