Hello,
I am performing a simulation in which I begin with FENE bonds and then switch to quartic bonds. I have read the quartic bonds page in the lammps manual, and understand the modified calculation for determining the energy after quartic bond failure. However, I am having an issue with the pressure as soon as I switch my system over to quartic (prior to breaking quartic bonds). I am wondering if this computational “sleight-of-hand” is also performed for the pressure tensor, because there seems to be an issue when bonds are changed from FENE to quartic–the simulation box quickly blows up. There are 2 log files written out in my simulation–1 for time lapsed under a FENE potential, and another for time lapsed for the quartic potential. These 2 log files contain a common timestep, which is the last timestep of FENE and the first timestep of quartic. Thus, their pressure values at this timestep should be equal, but the magnitude of the pressure in the quartic log file is significantly higher and the box rapidly expands at all times thereafter.
Here is the input file:
variable T equal 0.6
variable ts equal .01
variable EXPONENT equal 1/6
variable BASE equal (2^${EXPONENT})
units lj
dimension 3
atom_style bond
read_data …/lammps_inputs/interlocking_middle_beads_data_fix1.polymer
bond_style fene
bond_coeff 1 30.0 1.3 1.0 1.0
bond_coeff 2 30.0 1.3 1.0 1.0
*pair_style lj/cut {BASE}*
<i>pair_coeff * * 1.0 1.0 {BASE}
pair_coeff 3 4 3.0 1.0 ${BASE}
special_bonds lj/coul 0.0 1.0 1.0
velocity all create T {vseed} mom yes rot yes dist gaussian
fix 1 all press/berendsen iso 0.0 0.0 10
fix 2 all langevin $T T 2 {vseed}
fix 3 all nve
thermo_style custom step temp press pe ke ebond etotal epair lx ly lz xlo xhi ylo yhi zlo zhi
run_style respa 2 4 bond 1 pair 2
dump 1 all xyz 10000 …/generation/trajectory/trajectory_{usic}_{method}4_{k}_{trial}.gen.xyz
thermo 100
run 10000000
write_restart …/generation/restart/restart_{usic}_{method}4_{k}_{trial}.gen
bond_style quartic
bond_coeff 1 1200 -0.55 0.25 1.3 34.6878
bond_coeff 2 ${k} -0.55 0.25 1.3 34.6878
*pair_style lj/cut {BASE}*
<i>pair_coeff * * 1.0 1.0 {BASE}
pair_coeff 3 4 3.0 1.0 ${BASE}
special_bonds lj/coul 1.0 1.0 1.0
thermo_style custom step temp press pe ke ebond etotal epair lx ly lz xlo xhi ylo yhi zlo zhi
run_style respa 2 4 bond 1 pair 2
dump 2 all xyz 10000 …/quartic/trajectory/trajectory_{usic}_{method}4_{k}_{trial}.quar.xyz
thermo 100
run 1000000
write_restart …/quartic/restart/restart_{usic}_{method}4_{k}_{trial}.quar
Thank you for your time and advice.