Dear lammps users:
I am simulating a simple shear to a cubic bcc Fe. To do so, I fixed two end groups in the simulation box, and let the upper part have a small displacement gradually. The system is minimized and the stress is calculated at each step. I can get the correct strain-stress curve from the output of pxy. However, when I tried to calculate shear stress through force, Fx/Axz, the results become nonsense. I don’t know the reason, but I think the problem is that how to calculate the force. I have tried different methods like this:
boundary p s p
velocity lower set 0.0 0.0 0.0 units box
velocity upper set 0.0 0.0 0.0 units box
fix 10 lower setforce 0.0 0.0 0.0
fix 20 upper setforce 0.0 0.0 0.0
variable f atom fx
compute c1 mobile reduce sum v_f
compute c2 mobile group/group upper
variable fmx equal fcm(mobile,x)
variable stress1 equal (c_c1/lx/lz)
variable stress2 equal (c_c2[1]/lx/lz)
variable stress3 equal (v_fmx/lx/lz)
All these three stress values are different from the correct pxy output. Can you help me out? Thanks a lot!
Shijun Zhao
Greetings,
We can’t provide information on the error without your entire input script in this case since there are group declarations and fixes being omitted from the description.
Adrian
Thank you for the prompt reply. The whole input is the following if it helps:
units metal
dimension 3
boundary p s p
atom_style atomic
atom_modify map array
variable latparam1 equal 2.8665
lattice bcc ${latparam1}
region whole block 0 10 0 10 0 10 units lattice
create_box 1 whole
create_atoms 1 region whole
pair_style eam/fs
pair_coeff * * Fe.eam.fs Fe
variable ymax equal yhi
variable yup equal {ymax}-5.18
variable ymin equal ylo
variable ydown equal {ymin}+5.18
region lower block INF INF INF {ydown} INF INF units box
region upper block INF INF {yup} INF INF INF units box
group lower region lower
group upper region upper
group boundary union lower upper
group mobile subtract all boundary
min_style cg
minimize 1e-10 1e-10 5000 5000
variable f atom fx
compute c1 mobile reduce sum v_f
compute c2 mobile group/group upper
variable stress1 equal (c_c1/lx/lz)(1.61e-19/1e-30)1e-6
variable stress2 equal (c_c2[1]/lx/lz)(1.61e-19/1e-30)1e-6
variable fmx equal fcm(mobile,x)
variable stress3 equal (v_fmx/lx/lz)(1.61e-19/1e-30)*1e-6
thermo 100
thermo_style custom step temp etotal press pxx pyy pzz pxy pyz pxz v_stress1 v_stress2 v_stress3
velocity lower set 0.0 0.0 0.0 units box
velocity upper set 0.0 0.0 0.0 units box
fix 10 lower setforce 0.0 0.0 0.0
fix 20 upper setforce 0.0 0.0 0.0
variable pxy equal -pxy*0.1
label loop
variable a loop 100
displace_atoms upper move 0.01 0 0 units box
min_style cg
minimize 1e-10 1e-10 5000 5000
variable sx equal a*0.01/ly
print "XX {sx} {pxy} {stress1} {stress2} {stress3}"
next a
jump SELF loop
Shijun