Shearing a network

Dear all,
I need your help in using LAMMPS for bonded networks.
I am trying to reproduce the paper: Link.

I have generated a network with only bonded interactions. As per the paper, I want to apply a simple shear to the network and then do the second part of the simulation. So, when I applied a simple shear using change_box command and minimized it with the cg method, I ended up with a state with almost zero stress. I found the value of stress ~ 10^(-15).

LAMMPS version: LAMMPS (2 Aug 2023 - Update 3)
The input script:

#Set up basic simulation stuff
#=======================================================================================================================
dimension	2
units		lj
atom_style	molecular
atom_modify     id yes map hash sort 1 2.0
comm_modify	vel yes cutoff 2.0
boundary        p p p
neighbor	4.0 bin
neigh_modify	every 1 delay 0 check yes
read_data	initconf-Z4.1.d
mass		* 1
timestep	0.001

#Variables
#=======================================================================================================================
variable	L0		equal	$(lx)
variable	delta		equal	0.05*${L0}
variable	shear		equal	xy/${L0}
variable	sigma		equal	-pxy
velocity	all		create	0.0 2349851
reset_timestep	0

#Interaction parameters
#=======================================================================================================================
pair_style zero 4.0
pair_coeff * *

bond_style	harmonic/restrain
bond_coeff	1 1
run		0	

#Compute initial state
#=======================================================================================================================
dump		initdc  all custom 1 dump-Min-Init.lammpstrj id mol type x y vx vy fx fy 
dump_modify	initdc  format line '%d %d %d %0.16g %0.16g %0.16g %0.16g %0.16g  %0.16g' 
run		0
undump		initdc

print           "# step temp v_shear v_sigma" file Thermo.dat screen no

#Deform system [Simple shear]
#=======================================================================================================================
change_box	all triclinic
change_box	all xy delta ${delta} remap boundary pp pp pp units box
compute		gpres all pressure NULL bond virial # global pressure tensor (virial part)
run		0

dump		dc  all custom 1 dump-Min.lammpstrj id mol type x y vx vy fx fy 
thermo_style	custom step temp v_shear v_sigma pxx pyy pxy fmax
thermo_modify	press gpres
thermo		1

fix		1 all enforce2d
fix		2 all viscous 1
fix		3 all nve

min_style	cg
min_modify	dmax 1e-2 line forcezero #backtrack
minimize	1.0e-25 1.0e-12 100000000 100000000


print           "$(step) $(temp) $(v_shear) $(v_sigma)" append Thermo.dat screen no

undump		dc
uncompute	gpres
unfix		1 
unfix		2
unfix		3

2024-04-05T23:00:00Z

There is no question here.

Thanks for your reply akohlmey.

Apologies, I am sorry, I was not clear enough.

My question is, even after shearing the network [only bonded interaction and not allowed to break], the stress value is coming to be zero and I am expecting a non-zero number.
In the following, I am also attaching the data file for the system configuration. I want to ask is my minimization is correct.

initconf-Z4.1.d (7.3 KB)

Please note that LAMMPS doesn’t know anything about your system except for the geometry and potential choice and parameters in input. It will then compute forces and follow the forces to the (local) minimum geometry. As such the minimizer algorithm is correct since it has been tested and confirmed many, many times. However, what is open for discussion is whether the model you are using is different from your expectations. That is something that you have to work out for yourself or discuss with your adviser or supervisor or tutor or similar. People in the forum here don’t know about your research and what would be a suitable model, so cannot make an assessment of whether your input is suitable.

Thank you for your reply. I will follow your advice and discuss it locally with my group members.

Once again, thank you for your input comments. :slight_smile: