I would like to know if my stress output method is correct

When calculating a copper based nickel inclusion model at LAMMPS, I found that the output stress was not higher than that of the pure copper model, which should be abnormal. I am using the velocity command to shear the model, and I believe it may be due to a problem with the stress output code. Therefore, I am seeking help from my stress output code to see if it is correct.
Additionally, my model is developed through the one described in the “Dislocation and Inclusion” tutorial on the Atomsk official website. Specifically, I first create a hole in the copper substrate, then remove all unnecessary nickel atoms, leaving only a nickel sphere of similar size to the hole. Next, I move the nickel sphere into the hole and merge the two structures.

3d metal shear simulation

units metal
boundary s p p

atom_style atomic

read_data CuNi.lmp

pair_style eam/alloy
pair_coeff * * CuNi_v2.eam.alloy Cu Ni

min_style cg
minimize 1.0e-15 1.0e-15 10000 10000

neighbor 2.0 bin
neigh_modify delay 5

region lower block INF INF INF 20 INF INF units box
region upper block INF INF 180 INF INF INF units box
group lower region lower
group upper region upper
group boundary union lower upper
group mobile subtract all boundary

temp controllers

compute new3d mobile temp
compute new2d mobile temp/partial 0 1 1

compute myStress mobile stress/atom NULL
compute mobileStress mobile reduce sum c_myStress[4]
compute voro mobile voronoi/atom
compute v_total mobile reduce sum c_voro[1]

输出

variable step equal “step”
variable temp equal “temp”
variable press equal “press”
variable lx equal “lx”
variable l0 equal ${lx}
variable strain equal “(v_lx - v_l0)/v_l0”

variable tiji equal vol_mobile_region

variable v_mobile equal “group mobile sum c_v_atoms[1]”

variable v_mobile equal c_v_total
variable stress equal “c_mobileStress/c_v_total/10000”
thermo 100
fix def3 all print 500 “{step} {temp} {press} {strain} {stress} {v_mobile}” file data1.txt screen no

equilibrate

velocity mobile create 300.0 5812775 temp new3d

fix 1 all nvt temp 300.0 300.0 100
fix 2 boundary setforce 0.0 0.0 0.0

fix 3 mobile temp/rescale 10 300.0 300.0 10.0 1.0
fix_modify 3 temp new3d

thermo 100
thermo_modify temp new3d

timestep 0.001
run 100000
reset_timestep 0

shear

dump 1 all atom 1000 dump.shear

unfix 1

label next_step

velocity upper set 0.1 0 0 units box
velocity mobile ramp vx 0.0 0.1 y 20 180 sum yes units box

fix 3 mobile temp/rescale 10 300.0 300.0 10.0 1.0
fix_modify 3 temp new2d

thermo 100
thermo_modify temp new2d

run 100000

velocity upper set -0.1 0 0 units box
velocity mobile ramp vx 0 -0.1 y 20 180 sum yes units box
unfix 3
fix 3 mobile temp/rescale 10 300.0 300.0 10.0 1.0
fix_modify 3 temp new2d

thermo 100
thermo_modify temp new2d

run 200000

velocity upper set 0.1 0 0 units box
velocity mobile ramp vx 0 0.1 y 20 180 sum yes units box
unfix 3
fix 3 mobile temp/rescale 10 300.0 300.0 10.0 1.0
fix_modify 3 temp new2d
thermo 100
thermo_modify temp new2d

run 100000

if “${step} < 4000000” then “jump SELF next_step”

No serious scientists will take any results seriously from simulations that use fix temp/rescale.
It has been well established for a long time how bad and unsuitable this thermostat algorithm is for production calculations.

It is (obviously) syntactically correct. That is likely all the help you can get. This forum is not an input validation service, and it is the job of each individual researcher to validate their results.

Whether the results are meaningful is impossible to determine from remote without investing significant effort, having complete access to all your files, and being familiar with your research and its goals. It is not likely that people volunteering their time like the ones responding in a public forum will do that for you.

Thank you for your answer. I will continue to check my code.